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Abstract The theoretical investigation of charge (and spin) transport at 
nanometer length scales requires the use of advanced and powerful tech- 
niques able to deal with the dynamical properties of the relevant physical 
systems, to explicitly include out-of-cquilibrium situations typical for elec- 
trical/heat transport as well as to take into account interaction effects in a 
systematic way. Equilibrium Green function techniques and their extension 
to non-equilibrium situations via the Keldysh formalism build one of the pil- 
lars of current state-of-the-art approaches to quantum transport which have 
been implemented in both model Hamiltonian formulations and first-principle 
methodologies. In this chapter we offer a tutorial overview of the applications 
of Green functions to deal with some fundamental aspects of charge trans- 
port at the nanoscale, mainly focusing on applications to model Hamiltonian 
formulations. 
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The natural limitations that are expected to arise by the further miniaturiza- 
tion attempts of semiconductor-based electronic devices have led in the past 
two decades to the emergence of the new field of molecular electronics, where 
electronic functions are going to be performed at the single-molecule level, 
see recent overview in Refs. [TJ [2j El EJ [5j E] . The original conception which 
lies at the bottom of this fascinating field can be traced back to the paper by 
Ari Aviram and Mark Ratner in 1974 [7], where a single-molecule rectifying 
diode was proposed. Obviously, one of the core issues at stake in molecu- 
lar electronics is to clarify the question whether single molecules (or more 
complex molecular aggregates) can support an electric current. To achieve 
this goal, extremely refined experimental techniques are required in order to 
probe the response of such a nano-object to external fields. The meanwhile 
paradigmatic situation is that of a single molecule contacted by two metallic 
electrodes between which a bias voltage is applied. 

Recent experiments 

Enormous progress has been achieved in the experimental realization of 
such nano-devices, we only mention the development of controllable single- 
molecule junctions [8]-[22| and scanning tunneling microscopy based tech- 
niques [13]- [H]. With their help, a plethora of interesting phenomena like 
rectification [18] , negative differential conductance [9] [35] , Coulomb block- 
ade [23[l0im[T3[T!l2l], Kondo effect PUCE!], vibrational effects [10l[25l[ll 
[Ill[T6l[3ll[32l[33[35l[36l[2T],and nanoscale memory effects O EH SOJ US [44] , 
among others, have been demonstrated. 

The traditional semiconductor nanoelectronics also remains in front of 
modern research, in particular due to recent experiments with small quantum 
dots, where cotunncling effects were observed [45[ [46l [47], as well as new 
rectification effects in double quantum dots, interpreted as spin blockade J48[ 
[49l [50l [51] . Note, that semiconductor experiments are very well controlled at 
present time, so they play an important role as a benchmark for the theory. 

Apart from single molecules, carbon nanotubes have also found extensive 
applications and have been the target of experimental and theoretical studies 
over the last years, see Ref. [52] for a very recent review. The expectations 
to realize electronics at the molecular scale also reached into the domain of 
bio-molecular systems, thus opening new perspectives for the field due to the 
specific self-recognition and self-assembling properties of biomolcculcs. For 
instance, DNA oligomers have been already used as templates in molecular 
electronic circuits [53] [54] [55] . Much less clear is, however, if bio-molecules, 
and more specifically short DNA oligomers could also act as wiring systems. 
Their electrical response properties are much harder to disclose and there 



Green function techniques in the treatment of quantum transport 



3 



is still much controversy about the factors that determine charge migration 
through such systems [Ml UHl ESI EOl EH ESI EH - 

Theoretical methods 

The theoretical treatment of transport at the nanoscale (see introduction 
in EH ESI E3 ESI ESI ED]) requires the combined use of different tech- 
niques which range from minimal model Hamiltonians, passing through semi- 
empirical methods up to full first-principle methodologies. We mention here 
some important contributions, while we have no possibility to cite all relevant 
papers. 

Model Hamiltonians can in a straightforward way select, out of the many 
variables that can control charge migration those which are thought to be 
the most relevant ones for a specific molecule-electrode set-up. They con- 
tain, however, in a sometimes not well-controlled way, many free parameters; 
hence, they can point at generic effects, but they must be complemented 
with other methodologies able to yield microscopic specific information. Semi- 
empirical methods can deal with rather large systems due to the use of special 
subsets of electronic states to construct molecular Hamiltonians as well as to 
the approximate treatment of interactions, but often have the drawback of 
not being transferable. Ab initio approaches, finally, can deal in a very precise 
manner with the electronic and atomic structure of the different constituents 
of a molecular junction (metallic electrodes, molecular wire, the interface) 
but it is not apriori evident that they can also be applied to strong non- 
equilibrium situations. 

From a more formal standpoint, there are roughly two main theoretical 
frameworks that can be used to study quantum transport in nanosystems at 
finite voltage: generalized master equations (GME) [TTl [72] and nonequilib- 
rium Green function (NGF) techniques [73J [74j EH [76j ES] . The former also 
lead to more simple rate equations in the case where (i) the electrode-system 
coupling can be considered as a weak perturbation, and (ii) off-diagonal el- 
ements of the reduced density matrix in the eigenstate representation (co- 
herences) can be neglected due to very short decoherence times. Both ap- 
proaches, the GME and NGF techniques, can yield formally exact expressions 
for many observables. For non-interacting systems, one can even solve ana- 
lytically many models. However, once interactions are introduced - and these 
are the most interesting cases containing a very rich physics - different ap- 
proximation schemes have to be introduced to make the problems tractable. 

In this chapter, we will review mainly the technique of non-equilibrium 
Green functions. This approach is able to deal with a very broad variety 
of physical problems related to quantum transport at the molecular scale. 
It can deal with strong non-equilibrium situations via an extension of the 
conventional GF formalism to the Schwinger-Keldysh contour [74] and it 
can also include interaction effects (electron-electron, electron- vibron, etc) 
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in a systematic way (diagrammatic perturbation theory, equation of motion 
techniques). Proposed first time for the mesoscopic structures in the early 
seventies by Caroli et al. [ZZl EH [79l [80] , this approach was formulated in an 
elegant way by Meir, Wingreen and Jauho [5TJ [82l [HH [66l [84], who derived 
exact expression for noncquilibrium current through an interacting nanosys- 
tem placed between large nonintcracting leads in terms of the noncquilibrium 
Green functions of the nanosystem. Still, the problem of calculation of these 
Green functions is not trivial. We consider some possible approaches in the 
case of electron-electron and electron-vibron interactions. Moreover, as we 
will show later on, it can reproduce results obtained within the master equa- 
tion approach in the weak coupling limit to the electrodes (Coulomb block- 
ade), but it can also go beyond this limit and cover intermediate coupling 
(Kondo effect) and strong coupling (Fabry-Perot) domains. It thus offer the 
possibility of dealing with different physical regimes in a unified way. 

Now we review briefly some results obtained recently in the main directions 
of modern research: general nanoscale quantum transport theory, atomistic 
transport theory and applications to particular single-molecule systems. 

General nanoscale quantum transport theory 

On the way to interpretation of modern experiments with single-molecule 
junctions and STM spectroscopy of single molecules on surfaces, two main 
theoretical problems are to be solved. The first is development of appropriate 
models based on ab initio formulation. The second is effective and scalable 
theory of quantum transport through multilevel interacting systems. We first 
consider the last problem, assuming that the model Hamiltonian is known. 
Quantum transport through noninteracting system can be considered using 
the famous Landauer-Biittiker method [HU UHl HU IHH1 IH3 E01 EH EH , 
which establishes the fundamental relation between the wave functions (scat- 
tering amplitudes) of a system and its conducting properties. The method 
can be applied to find the current through a noninteracting system or through 
an effectively noninteracting system, for example if the mean-field description 
is valid and the inelastic scattering is not essential. Such type of an electron 
transport is called coherent, because there is no phase-breaking and quan- 
tum interference is preserved during the electron motion across the system. In 
fact, coherence is initially assumed in many ab initio based transport meth- 
ods (DFT+NGF, and others), so that the Landauer-Biittiker method is now 
routinely applied to any basic transport calculation through nanosystcms 
and single molecules. Besides, it is directly applicable in many semiconduc- 
tor quantum dot systems with weak electron-electron interactions. Due to 
simplicity and generality of this method, it is now widely accepted and is in 
the base of our understanding of coherent transport. 

However, the peculiarity of single-molecule transport is just essential role of 
electron-electron and electron-vibron interactions, so that Landauer-Biittiker 
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method is not enough usually to describe essential physics even qualitatively. 
During last years many new methods were developed to describe transport 
at finite voltage, with focus on correlation and inelastic effects, in particular 
in the cases when Coulomb blockade, Kondo effect and vibronic effects take 
place. 

Vibrons (the localized phonons) are very important because molecules 
arc flexible. The theory of clcctron-vibron interaction has a long history, 
but many questions it implies are not answered up to now. While the iso- 
lated electron-vibron model can be solved exactly by the so-called polaron 
or Lang-Firsov transformation [551 [Ml E2], the coupling to the leads pro- 
duces a true many-body problem. The inelastic resonant tunneling of single 
electrons through the localized state coupled to phonons was considered in 
Refs. [981 ISSl GUI [Ml QUI [103]. There the exact solution in the single- 
particle approximation was derived, ignoring completely the Fermi sea in the 
leads. At strong electron-vibron interaction and weak couplings to the leads 
the satellites of the main resonant peak are formed in the spectral function. 

The essential progress in calculation of transport properties in the strong 
electron-vibron interaction limit has been made with the help of the master 
equation approach [TOU QUI QUI QUI QUI Q03 Q3H1 QII1 032] ■ This method, 
however, is valid only in the limit of very weak molecule-to-lead coupling and 
neglects all spectral effects, which are the most important at finite coupling 
to the leads. 

At strong coupling to the leads and the finite level width the master equa- 
tion approach can no longer be used, and we apply alternatively the nonequi- 
librium Green function technique which have been recently developed to treat 
vibronic effects in a pcrturbative or self-consistent way in the cases of weak 
and intermediate electron-vibron interaction [1131 11141 11151 11161 11171 11181 

The case of intermediate and strong electron-vibron interaction at inter- 
mediate coupling to the leads is the most interesting, but also the most diffi- 
cult. The existing approaches are mean-field [1311I132|I133| . or start from the 
exact solution for the isolated system and then treat tunneling as a pertur- 
bation [TMQ35l[T3i[T3ZlQ3iQ3SQiQ]. The fluctuations beyond mean-field 
approximations were considered in Refs. [1411 1142] 

In parallel, the related activity was in the field of single-electron shuttles 
and quantum shuttles QUI EH EH EH EH QUI QUI EQl QSU QS21 Q53] . 
Finally, based on the Bardecn's tunneling Hamiltonian method [15411155111561 
11571 I158j and Tersoff-Hamann approach [1591 1160j , the theory of inelastic 
electron tunneling spectroscopy (IETS) was developed [1611 11621 11131 11141 
E3QI611I63]. 

The recent review of the electron-vibron problem and its relation to the 
molecular transport see in Ref. [164] . 

Coulomb interaction is the other important ingredient of the models, 
describing single molecules. It is in the origin of such fundamental effects as 
Coulomb blockade and Kondo efect. The most convenient and simple enough 
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is Anderson-Hubbard model, combining the formulations of Anderson impu- 
rity model [165] and Hubbard many-body model |166[ I167[ I168j . To analyze 
such strongly correlated system several complementary methods can be used: 
master equation and perturbation in tunneling, equation-of-motion method, 
self-consistent Green functions, renormalization group and different numeri- 
cal methods. 

When the coupling to the leads is weak, electron-electron interaction re- 
sults in Coulomb blockade, the sequential tunneling is described by the mas- 
ter equation method [1691 EQl HID QI21 QISl EH H23 HIS] and small co- 
tunneling current in the blockaded regime can be calculated by the next- 
order perturbation theory [1771 11781 1179] . This theory was used success- 
fully to describe electron tunneling via discrete quantum states in quantum 
dots [180(1181 1 182. 183J . Recently there were several attempts to apply master 
equation to multi-level models of molecules, in particular describing benzene 
rings [TM [15511156] . 

To describe consistently cotunncling, level broadening and higher-order 
(in tunneling) processes, more sophisticated methods to calculate the re- 
duced density matrix were developed, based on the Liouville - von Neumann 
equation [1871 QUI HH3 QSQl HSJ HMl [1931 MM or real-time diagrammatic 
technique [1941 [1951 EH HH1 G2E1 EH HH EE] ■ Different approaches were 
reviewed recently in Ref . [202] . 

The equation-of-motion (EOM) method is one of the basic and powerful 
ways to find the Green functions of interacting quantum systems. In spite of 
its simplicity it gives the appropriate results for strongly correlated nanosys- 
tems, describing qualitatively and in some cases quantitatively such impor- 
tant transport phenomena as Coulomb blockade and Kondo effect in quantum 
dots. The results of the EOM method could be calibrated with other available 
calculations, such as the master equation approach in the case of weak cou- 
pling to the leads, and the perturbation theory in the case of strong coupling 
to the leads and weak electron-electron interaction. 

In the case of a single site junction with two (spin-up and spin-down) states 
and Coulomb interaction between these states (Anderson impurity model), 
the linear conductance properties have been successfully studied by means of 
the EOM approach in the cases related to Coulomb blockade [203, 204 and 
the Kondo effect |205j . Later the same method was applied to some two-site 
models [2 ^ [207 ] [208 1 [209], Multi-level systems were started to be considered 
only recently [2101 1211j . Besides, there are some difficulties in building the 
lesser GF in the noncquilibrium case (at finite bias voltages) by means of the 
EOM method [2T^[2T3l[2T4l . 

The diagrammatic method was also used to analyze the Anderson impu- 
rity model. First of all, the perturbation theory can be used to describe weak 
electron-electron interaction and even some features of the Kondo effect [21 5j . 
The family of nonperturbative current-conserving self-consistent approxima- 
tions for Green functions has a long history and goes back to the Schwinger 
functional derivative technique, Kadanoff-Baym approximations and Hcdin 
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equations in the equilibrium many-body theory [2TB1 l21"7l 12151 I2T9I [2201 I22TI 
I222[ I223j . Recently GW approximation was investigated together with other 
methods l 224, 225, 226, 227 . It was shown that dynamical correlation effects 
and self-consistency can be very important at finite bias. 

Finally, we want to mention briefly three important fields of research, that 
we do not consider in the present review: the theory of Kondo effect [2281 
12291 [2301 12051 I23T1 [2321 [2331 [234] , spin-dependent transport [2351 EES [2371 
[2381 E39], and time-dependent transport [8% 1 [240 1 1241 ] l24"2 l [243] . 

Atomistic transport theory 

Atomistic transport theory utilizes semi-empirical (tight-binding [244) 1245] ) 
or ab initio based methods. In all cases the microscopic structure is taken 
into account with different level of accuracy. 

The most popular is the approach combining density-functional theory 
(DFT) and NGF and known as DFT+NGF [2461 EMU EH CSS [2501 EM 
[2521 I2%3l [2541 [2551 [2561 [2571 [2581 [2591 [2601 I26T1 [2621 [2631 [2641 [2651 [2661 
I267[ 1268 . This method, however is not free from internal problems. First 
of all, it is essentially mean-filed method neglecting strong local correlations 
and inelastic scattering. Second, density-functional theory is a ground state 
theory and e.g. the transmission calculated using static DFT eigenvalues will 
display peaks at the Kohn-Sham excitation energies, which in general do not 
coincide with the true excitation energies. Extensions to include excited states 
as in time-dependent density-functional theory, though very promising [269[ 
I270U271] , are not fully developed up to date. 

To improve DFT-based models several approaches were suggested, includ- 
ing inelastic electron-vibron interaction [RTl [2721 [273 QUE [27! HISl [2761 
[2771 [27811279 or Coulomb interaction beyond mean-field level [280] , or based 
on the LDA+U approachc [281] . The principally different alternative to DFT 
is to use an initio quantum chemistry based many-body quantum transport 
approach [Ml ESSE Ell EUS] . 

Finally, transport in bio-molecules attracted more attention, in particular 
electrical conductance of DNA [2861 EHIl EH3 [2891 [290], 

Outline 

The review is organized as follows. In Sections 2 we will first introduce 
the Green functions for non-interacting systems, and present few examples 
of transport through non-interacting regions. Then we review the master 
equation approach and its application to describe Coulomb blockade and 
vibron-mediated Franck-Condon blockade. In Section 3 the Keldysh NGF 
technique is developed in detail. In equilibrium situations or within the lin- 
ear response regime, dynamic response and static correlation functions are 
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related via the fluctuation-dissipation theorem. Thus, solving Dyson equa- 
tion for the retarded GF is enough to obtain the correlation functions. In 
strong out-of-equilibrium situations, however, dynamic response and corre- 
lation functions have to be calculated simultaneously and are not related 
by fluctuation-dissipation theorems. The Kadanof-Baym-Keldysh approach 
yield a compact, powerful formulation to derive Dyson and kinetic equations 
for non-equilibrium systems. In Sec. 4 we present different applications of the 
Green function techniques. We show how Coulomb blockade can be described 
within the Anderson-Hubbard model, once an appropriate truncation of the 
equation of motion hierarchy is performed (Sec. 4. A). Further, the paradig- 
matic case of transport through a single electronic level coupled to a local 
vibrational mode is discussed in detail within the context of the self-consistent 
Born approximation. It is shown that already this simple model can display 
non-trivial physics (Sec. 4.B). Finally, the case of an electronic system inter- 
acting with a bosonic bath is discussed in Sec. 4.C where it is shown that 
the presence of an environment with a continuous spectrum can modified the 
low-energy analytical structure of the Green function and lead to dramatic 
changes in the electrical response of the system. We point at the relevance 
of this situation to discuss transport experiments in short DNA oligomers. 
We have not addressed the problem of the (equilibrium or non-equilibrium) 
Kondo effect, since this issue alone would require a chapter on its own due 
to the non-perturbative character of the processes leading to the formation 
of the Kondo resonance. 

In view of the broadness of the topic, the authors were forced to do a 
very subjective selection of the topics to be included in this review as well as 
of the most relevant literature. We thus apologize for the omission of many 
interesting studies which could not be dealt with in the restricted space at 
our disposal. We refer the interested reader to the other contributions in this 
book and the cited papers. 
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2 From coherent transport to sequential tunneling 
(basics) 

2.1 Coherent transport: single-particle Green functions 

Nano-scalc and molecular-scale systems are naturally described by the discrete- 
level models, for example eigenstates of quantum dots, molecular orbitals, or 
atomic orbitals. But the leads are very large (infinite) and have a continuous 
energy spectrum. To include the lead effects systematically, it is reasonable 
to start from the discrete-level representation for the whole system. It can 
be made by the tight-binding (TB) model, which was proposed to describe 
quantum systems in which the localized electronic states play an essential 
role, it is widely used as an alternative to the plane wave description of elec- 
trons in solids, and also as a method to calculate the electronic structure of 
molecules in quantum chemistry. 

A very effective method to describe scattering and transport is the Green 
function (GF) method. In the case of non-interacting systems and coherent 
transport single-particle GFs are used. In this section we consider the matrix 
Green function method for coherent transport through discrete-level systems. 



(i) Matrix (tight-binding) Hamiltonian 

The main idea of the method is to represent the wave function of a particle 
as a linear combination of some known localized states i/> Q (r, cr), where a 
denote the set of quantum numbers, and a is the spin index (for example, 
atomic orbitals, in this particular case the method is called LCAO - linear 
combination of atomic orbitals) 

^(0=53^(0, (!) 

a 

here and below we use £ = (r, a) to denote both spatial coordinates and spin. 

Using the Dirac notations \a) = ip a (£,) and assuming that ip a (£.) are or- 
thonormal functions (a\(3) = 5 a p we can write the single-particle matrix 
(tight-binding) Hamiltonian in the Hilbcrt space formed by ipa{£,) 

H = J2( e »+e^ a )\a)(a\ + ^t a/9 |a)(/3|, (2) 

a a/3 

the first term in this Hamiltonian describes the states with energies e a , ip a 
is the electrical potential, the second term should be included if the states 
| a) are not eigenstates of the Hamiltonoian. In the TB model t a fj is the 
hopping matrix element between states \a) and |/3), which is nonzero, as a 
rule, for nearest neighbor cites. The two-particle interaction is described by 
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the Hamiltonian 



H = 



V M \am{5\{ 1 \, 



(3) 



in the two-particle Hilbert space, and so on. 

The energies and hopping matrix elements in this Hamiltomian can be 
calculated, if the single-particle real-space Hamiltomian is known: 



(4) 



This approach was developed originally as an approximate method, if the 
wave functions of isolated atoms are taken as a basis wave functions "0q(^)? 
but also can be formulated exactly with the help of Wannier functions. Only 
in the last case the expansion ((!]) and the Hamiltonian are exact, but 
some extension to the arbitrary basis functions is possible. In principle, the 
TB model is reasonable only when local states can be orthogonalizcd. The 
method is useful to calculate the conductance of complex quantum systems 
in combination with ah initio methods. It is particular important to describe 
small molecules, when the atomic orbitals form the basis. 

In the mathematical sense, the TB model is a discrete (grid) version of 
the continuous Schrodinger equation, thus it is routinely used in numerical 
calculations. 

To solve the single-particle problem it is convenient to introduce a new rep- 
resentation, where the coefficients c a in the expansion {T]) are the components 
of a vector wave function (we assume here that all states a are numerated 
by integers) 

C2 



(5) 



\cnJ 



and the eigenstates W\ are to be found from the matrix Schrodinger equation 



H^ A = E x $x, 

with the matrix elements of the single-particle Hamiltonian 



H, 



a/3 



e a + e<p a , a = (3, 
t a p, a / j3. 



(6) 



(7) 



Now let us consider some typical systems, for which the matrix method is 
appropriate starting point. The simplest example is a single quantum dot, the 
basis is formed by the eigenstates, the corresponding Hamiltonian is diagonal 
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Fig. 1 A finite linear chain of single-level sites. 
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The next typical example is a linear chain of single-state sites with only 
nearest-neighbor couplings (Fig. [I} 



H 



/ei ( ••• \ 
t e 2 t ■■■ 



• • • t tN-l t 

V • • • t e N J 



(9) 



The method is applied as well to consider the semi-infinite leads. Although 
the matrices are formally infinitely-dimensional in this case, we shall show 
below, that the problem is reduced to the finite-dimensional problem for the 
quantum system of interest, and the semi-infinite leads can be integrated out. 

Finally, in the second quantized form the tight-binding Hamiltonian is 



H = ^(e a + eip a ) clfia + X! t ^P C a 



C3. 



(10) 



(ii) Matrix Green functions and contact self-energies 

The solution of single-particle quantum problems, formulated with the help 
of a matrix Hamiltonian, is possible along the usual line of finding the 
wave- functions on a lattice, solving the Schrodingcr equation ([6]). The other 
method, namely matrix Green functions, considered in this section, was found 
to be more convenient for transport calculations, especially when interactions 
are included. 

The retarded single-particle matrix Green function G R (e) is determined 
by the equation 

[(e + in)l - H] G R = I, (11) 
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Fig. 2 A quantum system coupled to the left and right leads. 



where rj is an infinitesimally small positive number ?; = + . 

For an isolated noninteracting system the Green function is simply ob- 
tained after the matrix inversion 

G R = [(e + ir 1 )I-H}- 1 . (12) 

Let us consider the trivial example of a two-level system with the Hamiltonian 



H 



ei t 
t e 2 



The retarded GF is easy found to be (e = e 
Gi?(e) = (e-ei)(e-e 2 )+* 2 



irj) 



e - e 2 t 
t e-e 



(13) 



(14) 



Now let us consider the case, when the system of interest is coupled to 
two contacts (Fig. [2]). We assume here that the contacts are also described by 
the tight-binding model and by the matrix GFs. Actually, the semi- infinite 
contacts should be described by the matrix of infinite dimension. We shall 
consider the semi-infinite contacts in the next section. 

Let us present the full Hamiltonian of the considered system in a following 
block form 



(15) 



H 



h£ n LS 

n LS n S n RS 



H 



H, 



where , H^, and H 1 ^. arc Hamiltonians of the left lead, the system, and 
the right lead separately. And the off-diagonal terms describe system-to-lead 
coupling. The Hamiltonian should be hcrmitian, so that 



H 



SL — 



H 



LS> 



H 



SR 



H 



RS- 



(16) 



The Eq. (fTTj) can be written as 
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G L G LS 
G SL G s G SR | = I, 
Grs Gr 

where we introduce the matrix E = (e+ir])l, and represent the matrix Green 
function in a convenient form, the notation of retarded function is omitted in 
intermediate formulas. Now our first goal is to find the system Green function 
Gs which defines all quantities of interest. From the matrix equation (|17p 

(E - H° ) G LS - UlsGs = 0, (18) 
-H[ S G LS + (E - H° s ) G S - H RS G RS = I, (19) 
UrsGs + (E - H° ) Grs = 0. (20) 

From the first and the third equations one has 

G LS = {E-Uiy 1 H LS Gs, (21) 
Grs = (E — H R ) 1 H^Gs, (22) 

and substituting it into the second equation we arrive at the equation 

(E — H<j - E) G s = I. (23) 

where we introduce the contact self-energy (which should be also called re- 
tarded, we omit the index in this chapter) 

S = H[ s (E - Uiy 1 U LS + H RS (E - H^) 1 H RS . (24) 

Finally, we found, that the retarded GF of a nanosystem coupled to the 
leads is determined by the expression 

Gf(e)= [(e + ^I-Hl-S]" 1 , (25) 

the effects of the leads are included through the self-energy. 

Here we should stress the important property of the self-energy |24|) , it is 
determined only by the coupling Hamiltonians and the retarded GFs of the 
isolated leads G° R = (E - H^) ^ (i = L, R) 



S i = Hl s (E-H?) 1 n i s = U\ s G i R n iS , (26) 

it means, that the contact self-energy is independent of the state of the 
nanosystem itself and describes completely the influence of the leads. Later 
we shall see that this property conserves also for interacting system, if the 
leads are noninteracting. 

Finally, we should note, that the Green functions considered in this section, 
are single-particle GFs, and can be used only for noninteracting systems. 
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Fig. 3 A quantum system coupled to a semi-infinite ID lead. 



(iii) Semi-infinite leads 

Let us consider now a nanosystem coupled to a semi-infinite lead (Fig. [3]). 
The direct matrix inversion can not be performed in this case. The spectrum 
of a semi-infinite system is continuous. We should transform the expression 
(|26|) into some other form. 

To proceed, we use the relation between the Green function and the eigen- 
functions S'a of a system, which are solutions of the Schrodinger equation © . 
Let us define \P\(a) = c\ in the eigenstate |A) in the sense of definition ([5]), 
then 

fl&w-E* k(o) * !W) 

A 



e + it] — E\' 



(27) 



where a is the TB state (site) index, A denotes the eigenstate, Ex is the energy 
of the eigenstate. The summation in this formula can be easy replaced by the 
integration in the case of a continuous spectrum. It is important to notice, 
that the eigenfunctions \P\{a) should be calculated for the separately taken 
semi-infinite lead, because the Green function of isolated lead is substituted 
into the contact self-energy. 

For example, for the semi-infinite ID chain of single-state sites (n, m = 
1,2,...) 

^_ T dk* k {n)n{m) 



H H (A 



2tt 



IT] 



E k 



(28) 



with the eigenfunctions $fc(n) = V2sin kn, Ek = £o + 2tcosfc. 

Let us consider a simple situation, when the nanosystem is coupled only to 
the end site of the ID lead (Fig.[3|). From (|26|) we obtain the matrix elements 
of the self-energy 

Z a p = V* a V ip G™, (29) 

where the matrix element V\ a describes the coupling between the end site of 
the lead (n = m = 1) and the state \a) of the nanosystem. 
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To make clear the main physical properties of the lead self-energy, let us 
analyze in detail the semi-infinite ID lead with the Green function (f28|) . The 
integral can be calculated analytically ([70], p. 213, [244j ) 

li? _ 1 (* sin 2 kdk exp(iK(e)) 



Gf x (e) = -/ — - = (30) 

7T J _ JI £ + IT] — 6q — It COS k t 

K(e) is determined from e = 6q + 2tcosK. Finally, we obtain the following 
expressions for the real and imaginary part of the self-energy 

\Vla\ 



t 



ReS aa = (k - Vk 2 - 1 [6{k - 1) - 9{-k - 1)] ) , (31) 



Imi: aQ = J^t^Y^6{l - (32) 



k= C —^. (33) 
2t v ' 

The real and imaginary parts of the self-energy, given by these expressions, 
are shown in Fig.|H There are several important general conclusion that we 
can make looking at the formulas and the curves. 

(a) The self-energy is a complex function, the real part describes the energy 
shift of the level, and the imaginary part describes broadening. The finite 
imaginary part appears as a result of the continuous spectrum in the leads. 
The broadening is described traditionally by the matrix 

r = i (S - S f ) , (34) 

called level-width function. 

(b) In the wide-band limit (t — > oo), at the energies e — eo "C t, it is possible 
to neglect the real part of the self-energy, and the only effect of the leads is 
level broadening. So that the self-energy of the left (right) lead is 

= (35) 



(iv) Transmission, conductance, current 

After all, we want again to calculate the current through the nanosystem. 
We assume, as before, that the contacts are equilibrium, and there is the 
voltage V applied between the left and right contacts. The calculation of the 
current in a general case is more convenient to perform using the full power 
of the noncquilibrium Green function method. Here we present a simplified 
approach, valid for nonintcracting systems only, following Paulsson [291] . 

Let us come back to the Schrddinger equation ([6]) in the matrix represen- 
tation, and write it in the following form 
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Fig. 4 (Color) Real and imaginary parts of the contact self-energy as a function of energy 
for a one-band one-dimensional lead. 










= E\* S 







HL H° HL I I % I = E | % | , (36) 
Hrs H° 

where ^s, and are vector wave functions of the left lead, the nanosys- 
tem, and the right lead correspondingly. 

Now we find the solution in the scattering form (which is difficult to call 
true scattering because we do not define explicitly the geometry of the leads). 
Namely, in the left lead W L = + where is the eigenstate of H° , 
and is considered as known initial wave. The "reflected" wave as well 
as the transmitted wave in the right lead \Pr, appear only as a result of the 
interaction between subsystems. The main trick is, that we find a retarded 
solution. 

Solving the equation (|36|) with these conditions, the solution is 



# L ={1 + Gf H LS Gf H[ s ) (37) 

= G°f H RS G§ Hl s W° L (38) 
*s = Gf Hl s $l (39) 



The physical sense of this expressions is quite transparent, they describe the 
quantum amplitudes of the scattering processes. Three functions Wl, <?s, and 
\Pfi are equivalent together to the scattering state in the Landauer-Biittiker 
theory. Note, that Gj here is the full GF of the nanosystem including the 
lead self-energies. 

Now the next step. We want to calculate the current. The partial (for some 
particular eigenstate ^° A ) current from the lead to the system is 
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ji=L,R = - (^ f H lS % - 4Hj s ^) . (40) 

To calculate the total current we should substitute the expressions for 
the wave functions ([37|) - ([39|l . and summarize all contributions |291j . As a 
result the Landauer formula is obtained. We present the calculation for the 
transmission function. First, after substitution of the wave functions we have 
for the partial current going through the system 

h = Jl = -3r = -j (& R H RS &s - <4 H L^«) = 

- | (G°t - G° ) K RS G« Hi X) = 

J (^HwG^flGf H^lPg) . (41) 
The full current of all possible left eigenstates is given by 

' = = (< A H iS G^r fl Gf Ht s! pg A ) / l(E\), (42) 

A A 

the distribution function describes the population of the left states, 

the distribution function of the right lead is absent here, because we consider 
only the current from the left to the right. 

The same current is given by the Landauer formula through the transmis- 
sion function T(E) 

e f°° — 

I=-j T{E)f L {E)dE. (43) 

If one compares these two expressions for the current, the transmission 
function at some energy is obtained as 

T{E) = 2nJ2s(E - E x ) (*f A H^G^gM^) 

A 

= 27r ££<^ - E ^ KaHls^) (^G^GfH^) 

A <5 

= Tr (T L G$T R G§ ) . (44) 

To evaluate the sum in brackets we used the eigenfunction expansion (|2T|) for 
the left contact. 

We obtained the new representation for the transmission formula, which 
is very convenient for numerical calculations 

T = Tr (ffl) = Tr (T L G A T R G R ) . (45) 
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Finally, one important remark, at finite voltage the diagonal energies in 
the Hamiltonians H° , H° s , and H R are shifted e a — ► e a +eip a . Consequently, 
the energy dependencies of the self-energies defined by (|26|) are also changed 
and the lead self-energies are voltage dependent. However, it is convenient 
to define the self-energies using the Hamiltonians at zero voltage, in that 
case the voltage dependence should be explicitly shown in the transmission 
formula 

T(E) = Tr [T L (E - ecp L )G R (e)T R (E - e m )G A (e)] , (46) 

where ip R and ifL are electrical potentials of the right and left leads. 

With known transmission function, the current I at finite voltage V can 
be calculated by the usual Landauer-Biitiker formulas (without spin degen- 
eration, otherwise it should be multiplied additionally by 2) 

/oo 
T(E)[f L (E)-f R (E)]dE, (47) 

where the equilibrium distribution functions of the contacts should be written 
with corresponding chemical potentials fii, and electrical potentials (pi 



h(E) = — ^ ^ , f R (E) = j— ^ . (48) 

T 



cxp( g -y^ ) +l' cxp 
The zero- voltage conductance G is 



dl 

G - dV 



v=o 

where f°(E) is the equilibriumfcrmi- function 



- df°(E) 



cxp 



T 
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2.2 Interacting nanosy stems and master equation 
method 

The single-particle matrix Green function methods, considered in the previ- 
ous section, can be applied only in the case of noninteracting electrons and 
without inelastic scattering. In the case of interacting systems, the other ap- 
proach, known as the method of tunneling (or transfer) Hamiltonian (TH) , 
plays an important role, and is widely used to describe tunneling in supercon- 
ductors, in ferromagnets, effects in small tunnel junctions such as Coulomb 
blockade (CB), etc. The main advantage of this method is that it is easely 
combined with powerful methods of many-body theory. Besides, it is very 
convenient even for noninteracting electrons, when the coupling between sub- 
systems is weak, and the tunneling process can be described by rather simple 
matrix elements. 



2.2.1 Tunneling and master equation 

(i) Tunneling (transfer) Hamiltonian 

The main idea is to represent the Hamiltonian of the system (we consider 
first a single contact between two subsystems) as a sum of three parts: "left" 
H L , "right" H R , and "tunneling" H T 

H = H L +H R + H t , (51) 

Hl and Hr determine "left" \Lk) and "right" \Rq) states 

HlMO = EkMS), ( 52 ) 
Htfl> q (t) = E 9 1> q (S), (53) 

below in this lecture we use the index k for left states and the index q for right 
states. Ht determines "transfer" between these states and is defined through 
matrix elements Vk q = (Lk\HT\Rq) ■ With these definitions the single-particle 
tunneling Hamiltonian is 

H = Y,E k \k}(k\ + Y,E q \q)(Q\ + Yl [y^\d){k\ + V q \\k){q\] . (54) 

k£L qeR kq 

The method of the tunneling Hamiltonian was introduced by Bardeen 
|154j . developed by Harrison |155j . and formulated in most familiar second 
quantized form by Cohen, Falicov, and Phillips [156 . In spite of many very 
successful applications of the TH method, it was many times criticized for it's 
phenomcnological character and incompleteness, beginning from the work of 
Prange |157j . However, in the same work Prange showed that the tunneling 
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Hamiltonian is well defined in the sense of the perturbation theory. These 
developments and discussions were summarized by Duke [158] . Note, that 
the formulation equivalent to the method of the tunneling Hamiltonian can 
be derived exactly from the tight-binding approach. 

Indeed, the tight-binding model assumes that the left and right states can 
be clearly separated, also when they are orthogonal. The difference with the 
continuous case is, that we restrict the Hilbert space introducing the tight- 
binding model, so that the solution is not exact in the sense of the continuous 
Schrodingcr equation. But, in fact, we only consider physically relevant states, 
neglecting high-energy states not participating in transport. 

Compare the tunneling Hamiltonian ([54]) and the tight-binding Hamilto- 
nian ([2]), divided into left and right parts 



H=Y1 *M«X0l+ E ^><7l+ £ [V Sa \S)(a\+V s * a \a)(S\]. (55) 



The first two terms are the Hamiltonians of the left and right parts, the 
third term describes the left-right (tunneling) coupling. The equivalent matrix 
representation of this Hamiltonian is 



The Hamiltonians (|54|) and (|55|) arc essentially the same, only the first one 
is written in the eigenstate basis |fc), \q), while the second in the tight-binding 
basis | a), \(3) of the left lead and |<5), I7) of the right lead. Now we want to 
transform the TB Hamiltonian (|55[) into the eigenstate representation. 

Canonical transformations from the tight-binding (atomic orbitals) repre- 
sentation to the eigenstate (molecular orbitals) representation play an im- 
portant role, and we consider it in detail. Assume, that we find two unitary 
matrices Sl and S^, such that the Hamiltonians of the left part H° L and of 
the right part can be diagonalized by the canonical transformations 



Q/3GL 



a£L, SGR 




(56) 



m 




R — 



(57) 
(58) 



The left and right eigenstates can be written as 




(59) 



a: 




(60) 



6 



and the first two free-particle terms of the Hamiltonian (|54j) arc reproduced. 
The tunneling terms are transformed as 
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Hla, — S L HlrSr, (61) 

HL^S^HL^, (62) 

or explicitly 

J2 V 5a \S)(a\=J2V q k\q)(k\, (63) 

a£L,S£R kq 

where 

V q k = ^2 Vs a S Lak S R s q - (64) 

The last expression solve the problem of transformation of the tight-binding 
matrix elements into tunneling matrix elements. 

For applications the tunneling Hamiltonian (|54| should be formulated in 
the second quantized form. We introduce creation and annihilation Schrddinger 
operators c^ fe , CLk, <4j_, cn q . Using the usual rules we obtain 

H = H L ({4;ck}) +H R {{c\;c q }) +H T ({ci.;c fe ; C t; Cg }) , (65) 

H = ^(e fe + etp L (t))c\c h + ^(e 9 + e<p B {t))c\c q + ^ [v^cjcfc + V* k c\c q . 

k q kq 

(66) 

It is assumed that left c k and right c q operators describe independent 
states and are anticommutative. For nonorthogonal states of the Hamilto- 
nian Hi + Hr it is not exactly so. But if we consider Hl and Hr as two 
independent Hamiltonians with independent Hilbert spaces we resolve this 
problem. Thus we again should consider (|66j) not as a true Hamiltonian, but 
as the formal expression describing the current between left and right states. 
In the weak coupling case the small corrections to the commutation relations 
are of the order of \ V q k\ and can be neglected. If the tight-binding formulation 
is possible, (f6"6"|) is exact within the framework of this formulation. In general 
the method of tunneling Hamiltonian can be considered as a phenomeno- 
logical microscopic approach, which was proved to give reasonable results in 
many cases, e.g. in description of tunneling between superconductors and 
Josephson effect. 

(ii) Tunneling current 

The current from the state k into the state q is given by the golden rule 

27re 

J k ^ q = er qk = —\ Vqk \ 2 f L (k) (1 - f R (q)) 5{E k - E q ), (67) 
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the probability (1 — fn(E q )) that the right state is unoccupied should be 
included, it is different from the scattering approach because left and right 
states are two independent states! 

Then we write the total current as the sum of all partial currents from 
left states to right states and vice versa (note that the terms fL{k)fn{q) are 
cancelled) 

3 = it ^ |y " fc|2 [m ~ f{q)] 5{Eq ~ Ek) - (68) 

kq 

For tunneling between two equilibrium leads distribution functions are 
simply Fermi-Dirac functions (|48|) and current can be finally written in the 
well known form (To do this one should multiply the integrand on 1 = J 8{E— 
E q )dE.) 

e f 00 

J=-J T(E, V) [f L (E) - f R (E)} dE, (69) 

with 



T(E, V) = (2ir) 2 \Vk q \ 2 S(E - E k - e^ L )S{E - E q - &p R ). (70) 

qk 

This expression is equivalent to the Landauer formula (|47|) , but the trans- 
mission function is related now to the tunneling matrix clement. 

Now let us calculate the tunneling current as the time derivative of the 
number of particles operator in the left lead Nl = X)fc c l c fc- Current from 
the left to right contact is 



J(t) 



dN L 
dt 



le 
J 



H T ,N L 



(71) 



where is the average over time-dependent Schrodinger state. Nl com- 

mute with both left and right Hamiltonians, but not with the tunneling 
Hamiltonian 



H T ,N L 



[(Vqkclck + V qk C q cCj C {, 



(72) 



k' k q 

using commutation relations 

Ckc[,Ck' - c\,c k 'Ck = c k c ] k ,c k , + c\,c k c k > = (c k c\, + 5 kk , - c k c\,)c k > = 5 kk >c k , 
we obtain 



kq 



Vqk \ C \ C q 



(73) 



Now we switch to the Heisenberg picture, and average over initial time- 
independent equilibrium state 
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(6(t))=Sp(p eg d(t)), p eq = ^"ZZ^ ■ ( 74 ) 
One obtains 

J W = j E [ v «« <4(*) c *(*)> - v ,l (4(*K(*))] • ( 75 ) 

feq 

It can be finally written as 

j(t) = ^im fx;^p* a (*) ) = t Rc (E^ G ^>*) 



kq I \ kq 



We define " left-right" density matrix or more generally lesser Green func- 
tion 

Gt q {tut2) = i{c\{t2)c k {t x )). 

Later we show that these expressions for the tunneling current give the 
same answer as was obtained above by the golden rule in the case of nonin- 
teracting leads. 



(iii) Sequential tunneling and the master equation 

Let us come back to our favorite problem - transport through a quantum sys- 
tem. There is one case (called sequential tunneling), when the simple formulas 
discussed above can be applied even in the case of resonant tunneling 

Assume that a nonintcracting nanosystcm is coupled weakly to a thermal 
bath (in addition to the leads). The effect of the thermal bath is to break 
phase coherence of the electron inside the system during some time t p / 1; 
called decoherence or phase-breaking time. r p h is an important time-scale in 
the theory, it should be compared with the so-called "tunneling time" - the 
characteristic time for the electron to go from the nanosystem to the lead, 
which can be estimated as an inverse level- width function r . So that the 
criteria of sequential tunneling is 

r Tph < 1. (76) 

The finite decoherence time is due to some inelastic scattering mechanism 
inside the system, but typically this time is shorter than the energy relaxation 
time r e , and the distribution function of electrons inside the system can be 
nonequilibrium (if the finite voltage is applied), this transport regime is well 
known in semiconductor superlatticcs and quantum-cascade structures. 

In the sequential tunneling regime the tunneling events between the left 
lead and the nanosystem and between the left lead and the nanosystem are 
independent and the current from the left (right) lead to the nanosystcm 
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is given by the golden rule expression ([68]) . Let us modify it to the case of 
tunneling from the lead to a single level \a) of a quantum system 

3 = ¥ £ IKtfe|2 [/(fc) - Pa] 5{Ea ~ Ek) > (77) 

k 

where we introduce the probability P a to find the electron in the state \a) 
with the energy E a . 



(iv) Rate equations for noninteracting systems 

Rate equation method is a simple approach base on the balance of incoming 
and outgoing currents. Assuming that the contacts are equilibrium we obtain 
for the left and right currents 

Ji=L(R) = er ia [f?(E a ) - P a ] , (78) 

where 

fc 

In the stationary state J = Jl = —Jr, and from this condition the level 
population P a is found to be 

p r La fl{E a ) + rR a f R {E a ) 

Pa - p; — p , (oUj 

I La + I Rat 

with the current 

J = e r FLa f R r a (f L (E a ) f R {E a )). (81) 

It is interesting to note that this expression is exactly the same, as one can 
obtain for the resonant tunneling through a single level without any scatter- 
ing. It should be not forgotten, however, that we did not take into account 
additional level broadening due to scattering. 



(v) Master equation for interacting systems 

Now let us formulate briefly a more general approach to transport through in- 
teracting nanosystems weakly coupled to the leads in the sequential tunneling 
regime, namely the master equation method. Assume, that the system can 
be in several states |A), which are the eigenstates of an isolated system and 
introduce the distribution function P\ - the probability to find the system 
in the state |A). Note, that these states are many-particle states, for example 
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for a two-level quantum dot the possible states are |A) = 1 00} , |10), 01|), and 
|11). The first state is empty dot, the second and the third with one electron, 
and the last one is the double occupied state. The other non-electronic de- 
grees of freedom can be introduce on the same ground in this approach. The 
only restriction is that some full set of eigcnstates should be used 



(82) 



The next step is to treat tunneling as a perturbation. Following this idea, 
the transition rates r xx from the state A' to the state A are calculated using 
the Fermi golden rule 



rf* = - 
h 



2 



f\H T \i) 5{E f ~E t ). (83) 



Then, the kinetic (master) equation can be written as 

^ = £ rAVp A<-E rA ' A ^, (84) 

A' A' 

where the first term describes tunneling transition into the state |A), and the 
second term - tunneling transition out of the state |A). 

In the stationary case the probabilities arc determined from 

£ r xx 'p y = r y *p x . (85) 



For nonintcracting electrons the transition rates are determined by the 
single-electron tunneling rates, and are nonzero only for the transitions be- 
tween the states with the number of electrons different by one. For example, 
transition from the state |A') with empty electron level a into the state |A) 
with filled state a is described by 

r „ Q=1 „ a =o = r Laf 0( Ea) + r Ra f R (E a ), (86) 

where 2~£ Q and -Trq, are left and right level- width functions ([751) . 

For interacting electrons the calculation is a little bit more complicated. 
One should establish the relation between many-particle eigenstates of the 
system and single-particle tunneling. To do this, let us note, that the states 
|/) and |z) in the golden rule formula (|83|) are actually the states of the whole 
system, including the leads. We denote the initial and final states as 

\i) = \ki,X') = \ki)\\'), (87) 
\f) = \k f ,X) = \k f )\X), (88) 
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where k is the occupation of the single-particle states in the lead. The pa- 
rameterization is possible, because we apply the perturbation theory, and 
isolated lead and nanosystem are independent. 

The important point is, that the leads are actually in the equilibrium 
mixed state, the single electron states are populated with probabilities, given 
by the Fermi-Dirac distribution function. Taking into account all possible 
single-electron tunneling processes, we obtain the incoming tunneling rate 

F ™ =f E fi( E i*°) K**> X \Ht\ ik, \')\ 2 5(Ey + E lka -E X ), (89) 

ika 

where we use the short- hand notations: \ik, A') is the state with occupied 
fc-state in the i— th lead, while \ik, A) is the state with unoccupied fc-state in 
the i— th lead, and all other states are assumed to be unchanged, E\ is the 
energy of the state A . 

To proceed, we introduce the following Hamiltonian, describing single elec- 
tron tunneling and charging of the nanosystem state 

H T =J2 [Vxy k c k X xx ' + V^, k clx x ' x ] , (90) 

fcAA' 

the Hubbard operators X xx = |A)(A'| describe transitions between eigen- 
states of the nanosystem. 

Substituting this Hamiltonian one obtains 

r XX ' = yE fi( E ik.) \V lk „\ 2 \V xyk \ 2 5{E X , + E lka - E X ). (91) 

ika 

In the important limiting case, when the matrix element V\y k is k- 
independent, the sum over k can be performed, and finally 

r ™ = E r * ( Ex - E ^\ Vxx> i 2 ft ( Ex - Exi )• ( 92 ) 

i=L,R 

Similarly, the outgoing rate is 

roui = E r *( E * - E J \ v ^'\ 2 i 1 - - E *)) ■ ( 93 ) 

i=L,R 

The current (from the left or right lead to the system) is 

J,= L At) = e E (r&Px - r^ t Px) ■ (94) 

AA' 

This system of equations solves the transport problem in the sequential 
tunneling regime. 



Green function techniques in the treatment of quantum transport 27 

2.2.2 Electron-electron interaction and Coulomb blockade 



(i) Anderson-Hubbard and constant-interaction models 

To take into account both discrete energy levels of a system and the electron- 
electron interaction, it is convenient to start from the general Hamiltonian 

H = eafjd^d/3 + - ^ Vapasd^d^d^ds. (95) 

The first term of this Hamiltonian is a free-particle discrete-level model (|10p 
with e a p including electrical potentials. And the second term describes all 
possible interactions between electrons and is equivalent to the real-space 
Hamiltonian 

He, = \J d tJ dt'$H€$Ht')V(i,?)$(?)#(£), (96) 
where ip(£) are field operators 

= Z>«(Oda. (97) 

a 

ip a (£,) are the basis single-particle functions, we remind, that spin quantum 
numbers are included in a, and spin indices are included in £ = r, a as 
variables. 

The matrix elements are defined as 



v a0 , yS = j dzj der a (or^')v(ta^(OMe)- m 

For pair Coulomb interaction V(|r|) the matrix elements are 
V a ^ s = J2 fdcf dr'r a (r,vWp(r',<r')V(\r-r'\)^(r,a)Mr',<r')- (99) 



Assume now, that the basis states \a) are the states with definite spin 
quantum number a a . It means, that only one spin component of the wave 
function, namely i/) a ((T a ) is nonzero, and ipaio'a) = 0. In this case the only 
nonzero matrix elements are those with a a = <r 7 and erg = as, they are 

V a0 , lS = J dv J dv'r a {r)^W)V{\v - r'\)^(r)Mr')- (100) 

In the case of delocalized basis states tp a (r), the main matrix elements are 
those with a = 7 and (3 = 5, because the wave functions of two different 
states with the same spin are orthogonal in real space and their contribution 
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is small. It is also true for the systems with localized wave functions ip a (r), 
when the overlap between two different states is weak. In these cases it is 
enough to replace the interacting part by the Anderson-Hubbard Hamilto- 
nian, describing only density-density interaction 

Hah = ^ X! U a/3n a np. (101) 

a#/3 

with the Hubbard interaction defined as 

U a0 = JdvJ dr'\Mr)\ 2 \Mr')\ 2 V(\r - r'|). (102) 

In the limit of a single-level quantum dot (which is, however, a two-level 
system because of spin degeneration) we get the Anderson impurity model 
(AIM) 

Haim = 2J e <rdld a + UnfUi. (103) 

<r=n 

The other important limit is the constant interaction model (CIM), which 
is valid when many levels interact with similar energies, so that approxi- 
mately, assuming U a p = U for any states a and (3 

£r lv-n — U f^~Y U (\^~2\ UN(N-l) 
h ah = - U af3 n a np w — I ^ n a I - — I ^ n a I = . 

\ a } \ a ) 

(104) 

where we used h 2 = h. 

Thus, the CIM reproduces the charging energy considered above, and the 
Hamiltonian of an isolated system is 

Hcim = ^pdidp + E(N). (105) 

a/3 

Note, that the equilibrium compensating charge density can be easily in- 
troduced into the AH Hamiltonian 

Hah = ^ X! Ua < 3 ~ ("< 3 " ™< 9 ) ■ ( 106 ) 

a#/3 



(ii) Coulomb blockade in quantum dots 

Here we want to consider the Coulomb blockade in intermediate-size quantum 
dots, where the typical energy level spacing Ae is not too small to neglect it 
completely, but the number of levels is large enough, so that one can use the 
constant-interaction model (|105| . which we write in the cigenstate basis as 
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a 

where the charging energy E(n) is determined in the same way as previously, 
for example by the expression (|104|1 . Note, that for quantum dots the usage 
of classical capacitance is not well established, although for large quantum 
dots it is possible. Instead, we shift the energy levels in the dot e a = e a +e<p a 
by the electrical potential 

<Pa = V G + V R + Tl a (V L -V R ), (108) 

where rj a are some coefficients, dependent on geometry. This method can be 
easily extended to include any self-consistent effects on the mean-field level by 
the help of the Poisson equation (instead of classical capacitances). Besides, if 
all rj a are the same, our approach reproduce again the the classical expression 

Ecim = ^ £q?1q + E ( n ^ + en( Pext- (109) 

a 

The addition energy now depends not only on the charge of the molecule, 
but also on the state \a), in which the electron is added 

AE+ a (n,n a = 0->n+l,n a = 1) = E(n + 1) - E(n) + e a , (110) 

we can assume in this case, that the single particle energies are additive to the 
charging energy, so that the full quantum eigenstate of the system is |n,n), 
where the set h = {n a } shows weather the particular single-particle state |a) 
is empty or occupied. Some arbitrary state h looks like 

h = {n a } = (n l5 712,713, 714,715, ■■•) = (l, 1,0, 1,0, ...). (Ill) 

Note, that the distribution h defines also n = ^ Q n a . It is convenient, how- 
ever, to keep notation n to remember about the charge state of a system, 
below we use both notations \n, h) and short one \h) as equivalent. 

The other important point is that the distribution function f n (a) in the 
charge state |n) is not assumed to be equilibrium, as previously (this con- 
dition is not specific to quantum dots with discrete energy levels, the dis- 
tribution function in metallic islands can also be noncquilibrium. However, 
in the parameter range, typical for classical Coulomb blockade, the tunnel- 
ing time is much smaller than the energy relaxation time, and quasiparticlc 
noncquilibrium effects are usually neglected). 

With these new assumptions, the theory of sequential tunneling is quite 
the same, as was considered in the previous section. The master equation is 
[!8TlfT80lfT72lH82] 
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dj ^ A = E (nv-Mn -l,h',t) + r^p(n + l,n',t))- 

n' 

E " + rp* n ) P( n > ^ + I M". ^ ^ > ( 112 ) 

h' 

where p(n,n,t) is now the probability to find the system in the state \n,n), 
r^r 1 is the transition rate from the state with n—1 electrons and single level 
occupation h' into the state with n electrons and single level occupation n. 
The sum is over all states n', which are different by one electron from the state 
h. The last term is included to describe possible inelastic processes inside the 
system and relaxation to the equilibrium function p eq (n, n). In principle, it is 
not necessary to introduce such type of dissipation in calculation, because the 
current is in any case finite. But the dissipation may be important in large 
systems and at finite temperatures. Besides, it is necessary to describe the 
limit of classical single-electron transport, where the distribution function of 
qausi-particles is assumed to be equilibrium. Below we shall not take into 
account this term, assuming that tunneling is more important. 

While all considered processes arc, in fact, single-particle tunneling pro- 
cesses, we arrive at 



= E {^r™- x p{n,np - 0,t) + 5 n ^ n+1 p{h,n = l,i)) - 

E (s n0l r;- ln + s n0O r; +ln )p(n,t), (113) 


where the sum is over single-particle states. The probability p(h, np = 0, t) 
is the probability of the state equivalent to h, but without the electron in 
the state (3. Consider, for example, the first term in the right part. Here the 
delta- function 6 n g\ shows, that this term should be taken into account only if 
the single-particle state (3 in the many-particle state h is occupied, r^ n ^ 1 is 
the probability of tunneling from the lead to this state, p(h, rip = 0,i) is the 
probability of the state n', from which the system can come into the state n. 

The transitions rates are defined by the same golden rule expressions, as 
before, but with explicitly shown single-particle state a 



2tt 



C +1 " = yI' 



2n J2\V^\ 2 fk5(AE+ a -E k ), (114) 



h 

k 
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r 2a = Y \\ n ~ 1,n « = 0|#Ti|n,n a = 1^ 5{Ei-E f ) = 

Y E l^l 2 (! - /*) S ( AE n-l a Ek), (H5) 
fe 

there is no occupation factors (1 — f a ), f a because this state is assumed to 
be empty in the sense of the master equation (|113|) . The energy of the state 
is now included into the addition energy. 
Using again the level-width function 

r i=L , Ra (E) = ^ ]T \V lk , a \H(E - E k ). (116) 
fc 

we obtain 

rs +ln = r La f L (AE+ a ) + r Ra f° R {AE+ a ), (in) 
rr 1 " - r La (l - f° L (AE+_ la )) + r Ra (1 - f° R {AE+_ la )) . (us) 

Finally, the current from the left or right contact to a system is 
J l= L, R =eEEp(") r '« (<W°(AE+J ~ <W(1 - /°(AE+J)) . (119) 

a n 

The sum over a takes into account all possible single particle tunneling events, 
the sum over states n summarize probabilities p(h) of these states. 

(iii) Linear conductance 

The linear conductance can be calculated analytically [1811 1172] . Here we 
present the final result: 

G = y £ £ r2+ r rL Peg( "' na = 1} [1 _ f{AE ^ ' (120) 

a n—1 

where P eq (n,n a = 1) is the joint probability that the quantum dot contains 
n electrons and the level a is occupied 

P eq (n,n a = 1) = ^2p eq {n)S \n - ^ np 6 nal , (121) 

V / 

and the equilibrium probability (distribution function) is determined by the 
Gibbs distribution in the grand canonical ensemble: 
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Fig. 5 Linear conductance of a QD as a function of the gate voltage at different temper- 
atures T = O.QlEc, T = 0.03E C , T = 0.05-E c , T = 0.1E C , T = 0.15E C (lower curve). 



A typical behaviour of the conductance as a function of the gate volt- 
age at different temperatures is shown in Fig. [5j In the resonant tunneling 
regime at low temperatures T <§C At the peak height is strongly temperature- 
dependent. It is changed by classical temperature dependence (constant 
height) at T > At. 

(iv) Transport at finite bias voltage 

At finite bias voltage we find new manifestations of the interplay between 
single-electron tunneling and resonant free-particle tunneling. 

Now, let us consider the current-voltage curve of the differential conduc- 
tance (Fig. [7|). First of all, Coulomb staircase is reproduced, which is more 
pronounced, than for metallic islands, because the density of states is limited 
by the available single-particle states and the current is saturated. Besides, 
small additional steps due to discrete energy levels appear. This character- 
istic behaviour is possible for large enough dots with At <C Ec- If the level 
spacing is of the oder of the charging energy At ~ Ec, the Coulomb block- 
ade steps and discrete-level steps look the same, but their statistics (position 
and height distribution) is determined by the details of the single-particle 
spectrum and interactions [182] . 

Finally, let us consider the contour plot of the differential conductance 
(Fig. [7J). Ii is essentially different from those for the metallic island. First, 
it is not symmetric in the gate voltage, because the energy spectrum is re- 




(122) 
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Fig. 6 Coulomb staircase. 




Fig. 7 Contour plot of the differential conductance. 



stricted from the bottom, and at negative bias all the levels are above the 
Fermi-level (the electron charge is negative, and a negative potential means 
a positive energy shift). Nevertheless, existing stability patterns are of the 
same origin and form the same structure. The qualitatively new feature is 
additional lines correspondent to the additional discrete-level steps in the 
voltage-current curves. In general, the current and conductance of quantum 
dots demonstrate all typical features of discrete-level systems: current steps, 
conductance peaks. Without Coulomb interaction the usual picture of reso- 
nant tunneling is reproduced. In the limit of dense energy spectrum At — > 
the sharp single-level steps are merged into the smooth Coulomb staircase. 
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2.2.3 Vibrons and Pranck-Condon blockade 

(i) Linear vibrons 

Vibrons are quantum local vibrations of nanosystcms (Fig. [5]), especially im- 
portant in flexible molecules. In the linear regime the small displacements 
of the system can be expressed as linear combinations of the coordinates of 
the normal modes x q , which are described by a set of independent linear 
oscillators with the Hamiltonian 

^ = E (Sr g + \ m ^ &2 ) ■ (123) 

The parameters m q are determined by the microscopic theory, and p q 
(p q = —ifo]^r in the ^-representation) is the momentum conjugated to x q , 

[x q ,P q ]_ = ih. 

Let us outline briefly a possible way to calculate the normal modes of 
a molecule, and the relation between the positions of individual atoms and 
collective variables. We assume, that the atomic configuration of a system is 
determined mainly by the elastic forces, which are insensitive to the transport 
electrons. The dynamics of this system is determined by the atomic Hamil- 
tonian 

* I «t = Y,S; l +w{{Tln}h (124) 

n 

where W ({R„}) is the elastic energy, which includes also the static external 
forces and can be calculated by some ah initio method. Now define new 
generalized variables qi with corresponding momentum pi (as the generalized 
coordinates not only atomic positions, but also any other convenient degrees 
of freedom can be considered, for example, molecular rotations, center-of- 




Fig. 8 (Color) A local molecular vibration. The empty circles show the equilibrium posi- 
tions of the atoms. The energies e a , and the overlap integral t a fj are perturbed. 
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mass motion, etc.) 

Hat = J2^ + W{{qi}) ' (125) 

i 

"masses" m, should be considered as some parameters. The equilibrium co- 
ordinates q® are defined from the energy minimum, the set of equations is 

9W ({<#}) , x 

y * JJ =0. 126 

% 

The equations for linear oscillations are obtained from the next order ex- 
pansion in the deviations Aq^ = q$ — 

_ p? ^d 2 W({q°}) 

This Hamiltonian describes a set of coupled oscillators. Finally, applying 
the canonical transformation from Aqi to new variables x q {q is now the index 
of independent modes) 

i 

we derive the Hamiltonian (|123[) together with the frequencies uj q of vibra- 
tional modes. 

It is useful to introduce the creation and annihilation operators 



V, , (130) 



9 y/2 \y Ti 9 y/m q <jj q h i 
in this representation the Hamiltonian of free vibrons is {% = 1 



H^=J2^W ( 131 ) 



(ii) Elcctron-vibron Hamiltonian 

A system without vibrons is described as before by a basis set of states |a) 
with energies e a and inter-state overlap integrals £ Q( g, the model Hamiltonian 
of a noninteracting system is 

^s 0) = ( e « + et P»(t)) did* + ^pdldp, (132) 
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where (V a) d a are creation and annihilation operators in the states |a), and 
<p a (t) is the (self-consistent) electrical potential (|108j) . The index a is used 
to mark single-electron states (atomic orbitals) including the spin degree of 
freedom. 

To establish the Hamiltonian describing the interaction of electrons with 
vibrons in nanosystems, we can start from the generalized Hamiltonian 

Hs = ? « + E ^ ^ (133) 

where the parameters arc some functions of the vibronic normal coordinates 
x q . Note that we consider now only the electronic states, which were excluded 
previously from the Hamiltonian (|124p . it is important to prevent double 
counting. 

Expanding to the first order near the equilibrium state we obtain 

4» = £ £ ^*Ad a + ££ ^*A<b, (134) 

a q 9 a#/3 1 9 

where e Q (0) and f Q| a(0) are unperturbed values of the energy and the overlap 
integral. In the quantum limit the normal coordinates should be treated as 
operators, and in the second-quantized representation the interaction Hamil- 
tonian is 

= £ £ A «/ 3 K + 4)<& d f>- ( 135 ) 

a/3 q 

This Hamiltonian is similar to the usual electron-phonon Hamiltonian, but 
the vibrations are like localized phonons and q is an index labeling them, 
not the wave-vector. We include both diagonal coupling, which describes a 
change of the electrostatic energy with the distance between atoms, and the 
off-diagonal coupling, which describes the dependence of the matrix elements 
tap over the distance between atoms. 
The full Hamiltonian 

H = H° S + H V +H L + H R + H T (136) 

is the sum of the noninteracting Hamiltonian Hg, the Hamiltonians of the 
leads Hril), the tunneling Hamiltonian Ht describing the system-to-lead 
coupling, the vibron Hamiltonian Hy including electron-vibron interaction 
and coupling of vibrations to the environment (describing dissipation of vi- 
brons). 

Vibrons and the electron-vibron coupling are described by the Hamiltonian 

(n = i) 

j3\ a q + a q 

)dldp + H env . (137) 

q af3 q 
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The first term represents free vibrons with the energy hu q . The second term 
is the electron-vibron interaction. The rest part H env describes dissipation of 
vibrons due to interaction with other degrees of freedom, we do not consider 
the details in this chapter. 

The Hamiltonians of the right (R) and left (L) leads read as usual 

H l=L (R) = ^(dka + eipi)cl krT Ci k(r , (138) 
(pi are the electrical potentials of the leads. Finally, the tunneling Hamiltonian 

H T = E (Vika^C^da + V* k ^ a dic lka ) (139) 

i—L : R ka,a 

describes the hopping between the leads and the molecule. A direct hopping 
between two leads is neglected. 

The simplest example of the considered model is a single-level model 
(Fig. [9]) with the Hamiltonian 

H = e Q <$d + u (Ja+ A (a + + a) + ^ [^fc c 4 W + V ik4k d + h - c - » ( 14 °) 

ik 

where the first and the second terms describe free electron state and free vi- 
bron, the third term is electron-vibron interaction, and the rest is the Hamil- 
tonian of the leads and tunneling coupling (i = L, R is the lead index). 

The other important case is a center-of-mass motion of molecules between 
the leads (Fig.[T0|). Here not the internal overlap integrals, but the coupling to 
the leads Vika,a {x) is fluctuating. This model is easily reduced to the general 
model (|137[) . if we consider additionaly two not flexible states in the left and 
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right leads (two atoms most close to a system), to which the central system 
is coupled (shown by the dotted circles). 

Tunneling Hamiltonian includes ^-dependent matrix elements, considered 
in linear approximation 

H T = H (v ik<rta (x)4 k(T d a + fc.C.) , (141) 

V L , R (x) = V e^ L « V (lT^j ■ (142) 

Consider now a single-level molecule (a = 0) and extend our system, 
including two additional states from the left (a = I) and right (a = r) sides 
of a molecule, which arc coupled to the central state through x-depcndcnt 
matrix elements, and to the leads in a usual way through Then the 

Hamiltonian is of linear electron- vibron type 

Hm+v = 22 (ta + etp a ) d^da + ti(d]d + h.c.) + t r (dld + h.c.)+ 

a— Z,0,r 

+ ujqcJcl + (a + a*) ^A rfJdo - \i{d]do + h.c.) + \ r (d}.do + h.c.U . 

(143) 

(hi) Local polaron and canonical transformation 

Now let us start to consider the situation, when the electron-vibron interac- 
tion is strong. For an isolated system with the Hamiltonian, including only 



E 




Fig. 10 (Color) A center-of-mass vibration. 



Green function techniques in the treatment of quantum transport 39 

diagonal terms, 

Hs+v = ^d f a d a + ^ oj q a\a q X ^ a i + a «)4<k> ( 144 ) 

a q a q 

the problem can be solved exactly. This solution, as well as the method of the 
solution (canonical transformation) , plays an important role in the theory of 
electron- vibron systems, and we consider it in detail. 

Let's start from the simplest case. The single-level electron-vibron model 
is described by the Hamiltonian 



H, 



s+v 



= iocftd + ojna)a + A (a? + a) d'd, 



(145) 



where the first and the second terms describe free electron state and free 
vibron, and the third term is the electron-vibron interaction. 

This Hamiltonian is diagonalized by the canonical transformation (called 
"Lang-Firsov" or "polaron") [95l EH E7] 



with 



H = S^HS, 



A 



(a f - a) dU 

LUq 



(146) 
(147) 



S = exp 

the Hamiltonian Q145P is transformed as 

H s+V = S^Hs+vS = e Q d)d + Lu a^a + A (a* + a) fid, (148) 

it has the same form as f|145[) with new operators, it is a trivial consequence 
of the general property 



1 (fihh- 



and new single-particle operators are 

a = S~ 1 aS — a — —d^d, 
at = S-^S = at - ±dU, 
d= S^dS = exp [-A ( a t _ a )l d, 

Jt = g-irfS = exp [A (at - a)] d*. 



Substituting these expressions into (|148[) we get finally 



H. 



s+v 



= (i n --)dU 



w 



u)qcl' a. 



(149) 

(150) 
(151) 
(152) 

(153) 
(154) 
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We see that the electron- vibron Hamiltonian (|145[) is equivalent to the free- 
particle Hamiltonian (|154p . This equivalence means that any quantum state 
obtained as a solution of the Hamiltonian (|154p is one-to-one equivalent 
to the state as a solution of the initial Hamiltonian (|145p . with the same 
matrix elements for any operator 



It follows immediately that the eigenstates of the free-particle Hamiltonian 
are 

\ip nm ) = |n = 0,l;m = 0,l,2,...) = (dt)"i-^|0), (158) 



and the eigcn-cnergies are 



f=s- 1 fs, 



(155) 

(156) 
(157) 



E(n, 77i ) = ( ?o ) n + u>om. (159) 

The eigenstates of the initial Hamiltonian (1145j) are 

IVw) = S$ nm ) = e -^( at ~ a K Vr^=H0), (160) 

VTTl! 

with the same quantum numbers (n, m) and the same energies (|159p . This 
representation of the eigenstates demonstrates clearly the collective nature 
of the excitations, but it is inconvenient for practical calculations. 

Now let us consider the polaron transformation (|146p ~ p47p applied to the 
tunneling Hamiltonian 

Ht= Y, J2( V ^ c lk* d +V^Ci ka ) (161) 

i— L,i?. ka 

The electron operators in the left and right leads c^ a are not changed by 
this operation, but the dot operators d a , S a are changed in accordance with 
(|152p and (|153p . So that transformed Hamiltonian is 

ffr=EE (v ika e-^~ a )c\ k J + V; ka e^- a )Sc lka ) . (162) 

i—L^R ka 

Now we see clear the problem: while the new dot Hamiltonian (|154p is 
very simple and exactly solvable, the new tunneling Hamiltonian p62p is 
complicated. Moreover, instead of one linear electron- vibron interaction term, 
the exponent in (|162p produces all powers of vibronic operators. Actually, we 
simply remove the complexity from one place to the other. This approach 
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works well, if the tunneling can be considered as a perturbation, we consider 
it in the next section. In the general case the problem is quite difficult, but in 
the single-particle approximation it can be solved exactly [98l [99l I100( 1101] . 

To conclude, after the canonical transformation we have two equivalent 
models: (1) the initial model (|145[) with the eigenstates (|160p ; and (2) the 
fictional free-particle model f|154[) with the eigenstates (|158[) . We shall call 
this second model polaron representation. The relation between the models is 
established by (|155j) - ()157|) . It is also clear from the Hamiltonian (|148p . that 
the operators d' , d, , and a describe the initial electrons and vibrons in the 
fictional model. 



(iv) Inelastic tunneling in the single-particle approximation 

In this section we consider a special case of a single particle transmission 
through an electron-vibron system. It means that we consider a system cou- 
pled to the leads, but without electrons in the leads. This can be considered 
equivalently as the limit of large electron level energy eo (far from the Fermi 
surface in the leads). 

The inelastic transmission matrix T(e', e) describes the probability that 
an electron with energy e, incident from one lead, is transmitted with the 
energy e' into a second lead. The transmission function can be defined as the 
total transmission probability 

T{e) = J T{e' 1 e)de l . (163) 

For a noninteracting single-level system the transmission matrix is 

where r(e) = -Tt(e) + -0?( e ) is the level-width function, and A(e) is the real 
part of the self-energy. 

We can do some general conclusions, based on the form of the tunneling 
Hamiltonian (|162|) . Expanding the exponent in the same way as before, we 
get 



#r = E E V ** c Ld 

i=L,R ka \ 



q ° + E am (( flt ) m + a " 



with the coefficients 

A \ m e -(V"o) 2 /2 



(165) 



(166) 
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This complex Hamiltonian has very clear interpretation, the tunneling of one 
electron from the right to the left lead is accompanied by the excitation of 
vibrons. The energy conservation implies that 



(167) 



so that the inelastic tunneling with emission or absorption of vibrons is pos- 
sible. 

The exact solution is possible in the wide-band limit. [98l [99l flOOi 1101] 
It is convenient to introduce the dimensionless electron-vibron coupling 
constant 

At zero temperature the solution is 
T(e', e) = r L r R e- 2 <> V 9 —Me - e' - rnw Q ) 

z — ^ m< 



(168) 



3=0 



j\(m - j)\ l\ e - e + guj - {j + l)w + i-T/2 



(169) 



the total transmission function T(e) is trivially obtain by integration over e'. 
The representative results are presented in Figs. [TT1 and fT2l 

At finite temperature the general expression is too cumbersome, and we 
present here only the expression for the total transmission function 




Fig. 11 Transmission function as a function of energy at different electron-vibron coupling: 
g = 0.1 (thin solid line), g = 1 (dashed line), and g = 3 (thick solid line), at r = 0.1. 
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Fig. 12 Transmission function as a function of energy at different coupling to the leads: 
r = 0.01 (thin solid line), _T = 0.1 (dashed line), and _T = 1 (thick solid line), at g = 3. 



where n u is the equilibrium number of vibrons. 
(v) Master equation 

When the system is weakly coupled to the leads, the polaron representa- 
tion (|154p . (|162j) is a convenient starting point. Here we consider how the 
sequential tunneling is modified by vibrons. 

The master equation for the probability p(n, to, i) to find the system in 
one of the polaron eigenstates (|158|) can be written as 



where the first term describes tunneling transition into the state \n,m), and 
the second term - tunneling transition out of the state \n, m ), I v [p] is the 
vibron scattering integral describing the relaxation to equilibrium. The tran- 
sition rates -T^^, should be found from the Hamiltonian (|162|) . 

Taking into account all possible single-electron tunneling processes, we 
obtain the incoming tunneling rate 




dp(n, to) 
dt 



]T r^Pin', to') - J2 OX", m) + I V [ P ], (171) 




il 
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fi( E ik°) |(**> !. m \Ht\ ik, 0, m')| 2 <5(£ m' + E ika - E lm ) 



2tt 



J2fi( E ^) \Vik 





-^(a^-a) 




(m 







8 (Ea m ' + Eikcr — #lm) 



^ A(-Eito — Eom') \M mm '\ fi(Ei m — Eo m >), (172) 



i=L,R 



where 



e"» 



*-(<»+-<») 



(173) 



is the Franck-Condon matrix element. We use usual short-hand notations: 
\ik, n, m) is the state with occupied fc-state in the i— th lead, n electrons, and 
m vibrons, while \ik,n,m) is the state with unoccupied fc-state in the i— th 
lead, E nm is the polaron energy (|159|) . 
Similarly, the outgoing rate is 

^mm' = ^(Elmi — E 0m ) \M mm /\ 2 (l — fi(E\ m i — Eom)) ■ (174) 

The current (from the left or right lead to the system) is 

Ji= L>R {t) = e E (0(0, "0 " ^wf (1, m')) . (175) 



The system of equations (|171[) - ()175|) solves the transport problem in the 
sequential tunneling regime. 



(v) Franck-Condon blockade 

Now let us consider some details of the tunneling at small and large values 
of the electro- vibron coupling parameter g = I — 



The matrix element (|173[) can be calculated analytically, it is symmetric 
in m — m! and for m < m! is 

_ " (-gYV^n.e-^g^' -"0/2 

n(m-/)!(Z + m'-m)! " 1 ' 

The lowest order elements are 

m/2 

A/ 0m -e- 9/2 ^7=, (177) 

M n = (1 - 3 )e-^ 2 , (178) 
M 12 = ^(l-^)e^ 2 ... (179) 
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0.8^ — , — — — n 




Fig. 13 Franck-Condon matrix elements Mo m for weak (g = 0.1, squares), intermediate 
(g = 1, triangles), and strong (g = 10, circles) electron-vibron interaction. Lines are the 
guides for eyes. 




Fig. 14 The inverse life-time (tT)^ 1 as a function of \/u!o at optimal electron level 
position eo = A 2 /2ojo for neutral state (thin solid line), and for the charged state (dashed 
line), and for the neutral state at other level position eo = A 2 /4ljq (thick solid line). 



The characteristic feature of these matrix elements is so-called Franck- 
Condon blockade [1071 I108| . illustrated in Fig. [T3] for the matrix element 
Mo m . From the picture, as well as from the analytical formulas, it is clear, 
that in the case of strong electron-vibron interaction the tunneling with small 
change of the vibron quantum number is suppressed exponentially, and only 
the tunneling through high-energy states is possible, which is also suppressed 
at low bias voltage and low temperature. Thus, the electron transport through 
a system (linear conductance) is very small. 

There are several interesting manifestations of the Franck-Condon block- 
ade. 
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The life-time of the state |n,m) is determined by the sum of the rates of 
all possible processes which change this state in the assumption that all other 
states are empty 

n'm' 

As an example, let us calculate the life-time of the neutral state |0,0), 
which has the energy higher than the charged ground state |1, 0). 

T 00 = X] ^m'O = r i( Elm ~ ^°o) \M mQ \ 2 fi (Elm - E 00 ). (181) 

n'm' m i=L,R 

In the wide-band limit we obtain the simple analytical expression 

„m / \2 \ 

T oo= r ^- 9 h. f0 {^--+-o m ). (182) 

The corresponding expression for the life-time of the charged state (which 
can be excited by thermal fluctuations) is 

m / \2 \ 

The result of the calculation is shown in Fig. [Ml it is clear seen that the 
tunneling from the state |0, 0) to the charged state and from the state |1,0) 
to the neutral state is exponentially suppressed in comparison with the bare 
tunneling rate r at large values of the electron-vibron interaction constant 
A. This polaron memory effect can be used to create nano-memory and nano- 
switches. At finite voltage the switching between two states is easy accessible 
through the excited vibron states. It can be used to switch between memory 
states [TT2] . 

The other direct manifestation of the Franck-Condon blockade, - suppres- 
sion of the linear conductance, was considered in Refs. |107| 1108] , 
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3 Nonequilibrium Green function theory of transport 



3.1 Standard transport model: a nanosystem between 
ideal leads 

First of all, wc formulate a standard discrete-level model to describe nanoscalc 
interacting quantum systems (quantum dot, system of quantum dots, molecule, 
below " nanosystem" , " central system" , or simply " system" ) coupled to free 
conduction electrons in the leads. We include the Coulomb interaction with 
the help of the Anderson-Hubbard Hamiltonan to be able to describe cor- 
relation effects, such as Coulomb blockade and Kondo effect, which could 
dominate at low temperatures. At high temperatures or weak interaction 
the self-consistent mean-field effects are well reproduced by the same model. 
Furthermore, electrons are coupled to vibrational modes, below we use the 
electron- vibron model introduced previously. 

(i) The model Hamiltonian 

The full Hamiltonian is the sum of the free system Hamiltonian the 
inter-system electron-electron interaction Hamiltonian He, the vibron Hamil- 
tonian Hy including the electron-vibron interaction and coupling of vibra- 
tions to the environment (dissipation of vibrons), the Hamiltonians of the 
leads -fffl(L), and the tunneling Hamiltonian Ht describing the system-to- 
lead coupling 

H = H S + H C + H V + H L + H R + H T . (184) 

An isolated noninteracting nanosystem is described as a set of discrete 
states | a) with energies e a and inter-orbital overlap integrals t a p by the fol- 
lowing model Hamiltonian: 

Hs ] = Yl ( £ « + e V«(*)) 4,d a + (185) 

where d^ a ,d a arc creation and annihilation operators in the states \a), and 
<p a (t) is the effective (self-consistent) electrical potential. The index a is used 
to mark single-electron states (e.g. atomic orbitals) including the spin degree 
of freedom. In the eigenstate (molecular orbital) representation the second 
term is absent and the Hamiltonian is diagonal. 

For molecular transport the parameters of a model are to be determined 
by ab initio methods or considered as semi-empirical. This is a compromise, 
which allows us to consider complex molecules with a relatively simple model. 

The Hamiltonians of the right (R) and left (L) leads are 
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H l=L(R) = ^2(e lk a + eLp l {t))c i lk<J c tkrJ , (186) 

fc(T 

<fi(t) are the electrical potentials of the leads, the index k is the wave vector, 
but can be considered as representing an other conserved quantum number, 
a is the spin index, but can be considered as a generalized channel number, 
describing e.g. different bands or subbands in semiconductors. Alternatively, 
the tight-binding model can be used also for the leads, then (|186[) should be 
considered as a result of the Fourier transformation. The leads are assumed 
to be noninteracting and equilibrium. 
The tunneling Hamiltonian 

H T = Y, E (VikaMk^c + V^J^CfJ) (187) 
i—L,R kcr, a 

describes the hopping between the leads and the system. The direct hopping 
between two leads is neglected (relatively weak molecule-to-lead coupling 
case). Note, that the direct hoping between equilibrium leads can be easy 
taken into account as an additional independent current channel. 

The Coulomb interaction inside a system is described by the Anderson- 
Hubbard Hamiltonian 

Hc = o E UctfKnp. (188) 

This Hamiltonian is used usually only for the short-range part of Coulomb 
interaction. The long-range interactions can be better introduced through 
the self-consistent electrical potential ip a , which is determined by the Poison 
equation with the average electron density. 

Vibrations and the electron-vibron coupling are described by the Hamil- 
tonian 

H v = ]T Rwgojo, A y°9 + a i)<& d f> + He- (189) 

q a/3 q 

Here vibrations are considered as localized phonons and q is the index labeling 
them, not the wave-vector. The first term describes free vibrons with the 
energy TiLo q . The second term represents the electron-vibron interaction. The 
third term describes the coupling to the environment and the dissipation of 
vibrons. We include both diagonal coupling, which originates from a change 
of the electrostatic energy with the distance between atoms, and the off- 
diagonal coupling, which can be obtained from the dependence of the matrix 
elements t a p over the distance between atoms. 

(ii) Nonequilibrium current and charge 

To connect the microscopic description of a system with the macroscopic 
(elcctrodynamic) equations and calculate the obscrvables, we need the expres- 
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sions for the nonequilibrium electrical charge of the system and the current 
between the system and the leads. 

The charge in a nonequilibrium state is given by (Qo is the background 
charge) 

Qs(t) = eY / (dU a )-Qo- (190) 

a 

To calculate the current we find the time evolution of the particle number 
operator Ns = ^2 a d' a d a due to tunneling from the left (i = L) or right 
(i = R) contact. 

The current from the left (i = L) or right (i = R) contact to the nanosys- 
tem is determined by (note, that we consider e as the charge of the electron 
(negative) or the hole (positive)) 

where 

H t = E (V ika<a 4 ka d a + V* k ^ a dlc lka ) (192) 

kcr, a 

is the Hamiltonian of the coupling to the corresponding contact. The current 
is determined by this only part of the full Hamiltonian (|136fl . because all 
other terms commute with 7Vs. 

Applying the commutation relation 

d n , d a d 



dadtdp — dj^dpda =d a d { pdfj + d^dadp = 
{dadfa + 6 a p - d a cfp)dp = 6 a pd a , (193) 

one obtains finally 

= J E [W ( c l*d a ) - V* k ^ a (dic lka )] . (194) 



a. a 



(iii) Density matrix and NGF 

The averages of the operators in Eqs. (|190|) and p94j) are the elements of the 
density matrix in the single-particle space 

P* a (t) = (4(t)d a (t)> , (195) 
p a ,ika(t) = (c\ k(T (t)d a (t)) . (196) 

It is possible, also, to express it as a two-time Green function at equal 
time's 
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,(t) =eJ2p aa (t) = -ie^G< a (t,t), (197) 



Ji(t) = ^Im Vika, a Pa,ikc(t) = ^Re V ^ a G< lka {t,t) , 

\k<r,a J \k<r,a J 

(198) 

where we define the system-to-lead lesser Green function 

G< ika (ti,t2) = i{4 ka (h)d a (ti)), (199) 

while noncquilibrium charge distribution of the molecule is determined by 
the system lesser function 

G< p (h,t2) = i (4(*a)da(*i)) • (200) 

One can ask: what is the advantage to use the more complex two-time 
Green functions instead of density matrices? There are several reasons. First 
of all, NGF give, as we shall see below, a clear description of both density 
of states and distribution of particles over this states. Then, the equations 
of motion including interactions and the influence of environment can be 
obtained with the help of a diagrammatic technique, and (very important) 
all diagrammatic results of equilibrium theory can be easily incorporated. 
Retardation effects are conveniently taken into account by two-time Green 
functions. And, ... finally, one can always go back to the density matrix when 
necessary. 

It is important to note, that the single-particle density matrix (|195p should 
not be mixed up with the density matrix in the basis of many-body eigenstates. 

In these review we consider different methods. The density matrix can be 
determined from the master equation. For Green functions the EOM method 
or Keldysh method can be applied. Traditionally, the density matrix is used 
in the case of very weak system-to-lead coupling, while the NGF methods 
are more successful in the description of strong and intermediate coupling to 
the leads. The convenience of one or other method is determined essentially 
by the type of interaction. Our aim is to combine the advantages of both 
methods. 
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3.2 Nonequilibrium Green functions: definition and 
properties 

In the previous section we found, that the current through a system (as well 
as other observables) can be expressed through nonequilibrium Green func- 
tions. Here we give the definitions of retarded, advanced, lesser, and greater 
Green functions and consider some simple examples. We also introduce a very 
important concept of the Schwinger-Keldysh closed-time contour, and define 
contour Green functions. This section is a little bit technical, but we need 
these definitions in the next sections. 



3.2.1 Spectral - retarded (G R ) and advanced (G A ) functions 

(i) Definition 

Retarded Green function for fermions is defined as 

G^(t 1 ,t 2 ) = -te(t 1 -t 2 )(\c Q (t 1 ),cl{t 2 )] ), (201) 



where c£, (t) , c a (t) are creation and annihilation time-dependent (Heisenberg) 
operators, [c, d]+ = cd + dc is the anti-commutator, (...) denotes averaging 
over equilibrium state. 

We use notations a, (3, ... to denote single-particle quantum states, the 
other possible notation is more convenient for bulk systems 



G R (x 1 ,x 2 ) = -i0(h -t 2 )([c( a;i ),c t (x 2 )] + ), 



(202) 



where x = r, t, a, ... or x = k, t, a, etc. Some other types of notations can 
be found in the literature, they are equivalent to (|20 1 1) . 
The advanced function for fermions is defined as 



(203) 



Finally, retarded and advanced functions for bosons can be defined as 



G%{ti,h) = -W(h - t 2 ) ^ [a a (ti),aj(f a ) 
G%{tiM = i6{t 2 - t t ) ( [a a (ti),aj(f 2 ) 



(204) 



(205) 



where (V a (t), a a (t) are creation and annihilation boson operators, [a, 6]_ 
ab — ba is the commutator. 
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(ii) Discussion of averaging 

The average value of any operator O can be written as (O) = (t\0 \t) in 
the Schrodinger representation or (O) = (0\O (t)\0) in the Heisenberg rep- 
resentation, where |0) is some initial state. This initial state is in principle 
arbitrary, but in many-particle problems it is convenient to take this state as 
an equilibrium state, consequently without time-dependent perturbation we 
obtain usual equilibrium Green functions. 

In accordance with this definition the Heisenberg operators c a (t), ct(t), 
etc. are equal to the time-independent Schrodinger operators at some initial 
time t$: c a (to) = c a , etc. Density matrix of the system is assumed to be 
equilibrium at this time p(to) = p e q- Usually we can take io = for simplicity, 
but if we want to use to the transformation to Heisenberg operators 
should be written as 

f H (t) = e iH(t-to)fs e -iH(t-to), (206 ) 

In fact, the initial conditions are not important because of dissipation 
(the memory about the initial state is completely lost after the relaxation 
time). However, in some pathological cases, for example for free noninter- 
acting particles, the initial state determines the state at all times. Note also, 
that the initial conditions can be more convenient formulated for Green func- 
tions itself, instead of corresponding initial conditions for operators or wave 
functions. 

Nevertheless, thermal averaging is widely used and we define it here explic- 
itly. If we introduce the basis of exact time-independent many-particle states 
|n) with energies E n , the averaging over equilibrium state can be written as 

<6> = i£ -- EJT (» O h (t) n), Z = ^ e- E » /T - (207) 



In the following when we use notations like (O) or ( \P 



0(t) 



if - ), we as- 



sume the averaging with density matrix (density operator) p 

[6) = Sp (p6) , (208) 



for equilibrium density matrix and Heisenberg operators it is equivalent to 
([2071). 



(iii) Free-particle retarded function for fermions 



Now consider the simplest possible example - retarded Green function for 
free particles (fermions). 
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The free-particle Hamiltonian has equivalent form if one uses Schrodinger 
or Heisenberg operators 



(209) 



because (here we assume t = 0) 



t f+\ fj.\ iHt t —iHt iHt -iHt 

c a (t)c a (t) = e c> a e e c a e 

iHt J „ -iHt _ J 



(210) 



where we used that c^Cq, is commutative with the Hamiltonian H = ^ Q e a <^ a c a . 
From the definitions (|20Tj) and ([207ll 



c a (ti),4(t 2 )] + \ = (c a (t 1 )cJ(t 2 )+4(t 2 )c Q (ti)J> = 
h)e- i6tl e i6t2 cl(t 2 )e- iAt2 + e ltit \ 



(211) 



G%p(h,h) = -i6(h - f a ) / [c a (fi), cJ(t 2 )] + \ 



(212) 



where we used some obvious properties of the creation and annihilation op- 
erators and commutation relations. 

We consider also the other method, based on the equations of motion for 
operators. From Liuville - von Neuman equation we find (all c-operators are 
Heisenberg operators in the formula below, (f) is omitted for shortness) 



i^^- = [c a (t),H]_=^2ei3 \c a ,c^cp 



dt 



13 13 

= ^£/3 (c a Cp + CpCa^J cp — epSgpcp = e a c a (f), 

13 13 

so that Heisenberg operators for free fermions are 

c Q (f) = e-"°* Cct (0), C t(i) = e*«* c t (0). 



(213) 



(214) 
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Substituting these expressions into (|201|) we obtain again (]212|) . Note also 
that if we take to ^ 0, then Heisenberg operators for free fermions are 

c a (t) = e-^'-^Cafo), 4 (t) = e «"('-*o) c t (f (J ) | (215) 
but the result for the Green functions is just the same, because 

Ca(ti),C^(i 2 ) + ^ = ^C a (ti)c^(t 2 ) + C^(t 2 )c a (ti)^ 

= e M*2-to)-«a(ti-to) / CaC t + c t ^ = e-^-C*!-*")^. (216) 

It is interesting to make Fourie-transform of this function. In equilibrium 
two-time function G R a(ti, t 2 ) is a function of the time difference only, so that 
wc define transform over time difference (t\ — t 2 ) 

/•OO 

G R (e) = / G R {h - t 2 )e 4(e+40)(tl ^ t2) d(ti - f 2 ), (217) 
Jo 

we add infinitely small positive complex part to e to make this integral well 
defined in the upper limit (this is necessary for free particles without dissi- 
pation because function (|2 1 2(1 oscillates at large times t = t\ — t 2 and the 
integral (2"T7jl can not be calculated without iO term. Then we obtain 



= — (218) 



More generally, transformation (|217[) can be considered as the Laplas 
transformation with complex argument z = e + irj. 
For advanced function 

G%{t 1 M)=iO{t 2 -t l )e- ie ^- t ^8 a p, (219) 

the Fourier transform is given by 

r° 

G A {e)= G A (t 1 -t 2 )e^ e - lQ) ^- t ^d(t 1 ~t 2 ), (220) 

J — oo 

with other sign of the term iO. 



(iv) Spectral function 

Finally, we introduce the important combination of retarded and advanced 
functions known as spectral or spectral weight function 

A aP (e)=i(G^(e)-G A (e)), (221) 
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in equilibrium case Fourie-transformed retarded and advanced functions are 
complex conjugate G A (e) = (G R (e))*, and A a p(e) = -21mG^Je). 
For free fermions the spectral function is 

A aP {e) = -21m ( ^7^) = 2vr5(e - e a )5 a0 . (222) 

\e - e a + lOJ 

The result is transparent - the function A a p(e) is nonzero only at particle 
eigcn-cnergies, so that 

p(e) = ^-SpA^(e) = i- ]T A aa {e) =^S(e- e a ) (223) 
Zir Zir * — ' 1 — ' 

a a 

is the usual energy density of states. Note that the imaginary part iO is 
necessary to obtain this result, thus it is not only mathematical trick, but 
reflects the physical sense of the retarded Green function. 
If we introduce finite relaxation time 

O r ) = -iB^e-^-^Sap, (224) 
then the spectral function has familiar Lorcntzian form 



A M = 7 21 2 - ( 22 5) 
(e - e a y + Y 



Finally, spectral function has a special property, so-called sum rule, namely 

A a p{e)^- = S af} . (226) 



3.2.2 Kinetic - lesser (G < ) and greater (G > ) functions 

(i) Definition 

Spectral functions, described before, determine single-particle properties of 
the system, such as quasiparticle energy, broadening of the levels (life-time), 
and density of states. These functions can be modified in noncquilibrium 
state, but most important kinetic properties, such as distribution function, 
charge, and current, are determined by lesser Green function 

G< /3 (t 1 ,t 2 )=i(cl(t 2 )c a (t 1 )). (227) 

Indeed, density matrix is the same as equal-time lesser function 

p a p(t) = (cj(t)ca(t)) = -iG< (t,t). (228) 
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the number of particles in state |a) (distribution function) is 

»«(*) = (cUt)c a (t)) = -iG< a (t,t), (229) 
the tunneling current is 

j w = *i E K (4(*)^(*)) - ^ (4(*K(*)) 



kq 



2e 



Re \J2v qk G<(t,t) 

kq 



In addition to the lesser the other (greater) function is used 

G^ J (ti ) t2) = -i(c a (ti)cJ(t 2 )). 
For bosons lesser and greater functions are defined as 

G< l3 (t 1 ,t 2 )=-i(a} j (t 2 )a a (t 1 ) 



G^{t 1 ,t 2 ) = -i(a ol {t l )a\{t 2 )) 



(230) 



(231) 



(232) 



(233) 



The name "lesser" originates from the time-ordered Green function, the 
main function in equilibrium theory, which can be calculated by diagrammatic 
technique 

G af3 (h,h) = -*(T(c a (ti)cJ(t a ))), (234) 



G a p(ti, t 2 ) 



-i (c a (ti)cp(t 2 )\ if ii>ta 
z^(i 2 )c a (ii)\ if ii<t 2 



G Q /3 = G > 
=> G Q /3 = G 



Q/3' 



(235) 



Q/3' 



here additional sing minus appears for interchanging of fcrmionic creation- 
annihilation operators. Lesser means that t\ < t 2 . 

From the definitions it is clear that the retarded function can be combined 
from lesser and greater functions 



Gaf3{t\,t 2 ) — 6(t\ - t 2 ) G>p(ti,t 2 ) - G<p(tl,fa) 



(236) 



(ii) Free-particle lesser function for fermions 

Now let us consider again free fermions. Heisenbcrg operators for free fermions 
are (t = 0) 

Ca (t) = e-"°*c Q (0), c t( t ) = e *«* c t (0). (237) 
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Lesser function is 

GipfaM =i(cl{t 2 )c a {h))=ie i ^- i ^ (c\c a ) 

= ie- ie ^-^f o (e a )S a0 , (238) 

one sees that contrary to the retarded function, the lesser function is propor- 
tional to the distribution function, in equilibrium this is Fermi distribution 
function 

= . (239) 

e t +1 

It is interesting to compare this answer with the result for nonthermal 
initial conditions. Assume that initial state is described by the density matrix 

Pa/3 = (^i3 c a^ '■> now w ith nonzero off-diagonal elements. Time dependence of 
the density matrix is given by 

Pap® = e^-^W (240) 

We obtain the well known result that off-diagonal elements oscillate in time. 
Now define Fourier-transform for lesser function (r = t\ —t<i) 

/>oc 

G<(e)= G < (r)e l [ e+i0sign ( r )] r dr, (241) 



note that here we use Fourie-transform with complicated term iOsign(r), 
which makes this transformation consistent with previously introduced trans- 
formations (|217p for retarded (r > 0) and (I220|) advanced (r < 0) functions. 
Applying this transformation to (|238[) we obtain 



J — OC 

= 2Trif°(e a )5(e~e a )6 a p. (242) 
For free fermion greater function one obtaines 

GlpihM) = -le-^-^il - f°(e a ))6 af3 , (243) 
G> (e) = -2%i(l - f(e a ))8(e - e a )5 a p. (244) 



(iii) Equilibrium case. Fluctuation-dissipation theorem 

Now we want to consider some general properties of interacting systems. In 
equilibrium the lesser function is not independent and is simply related to 
the spectral function by the relation 
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G<,(e) = iA a p{e)f{e). (245) 

This relation is important because establish equilibrium initial condition for 
noncquilibrium lesser function, and propose useful Ansatz if equilibrium dis- 
tribution function / (e) is replaced by some unknown noncquilibrium func- 
tion. 

Here we prove this relation using Lehmann representation - quite useful 
method in the theory of Green functions. The idea of the method is to use 
exact many-particle eigenstates |n), even if they are not explicitly known. 

Consider first the greater function. Using states \n) we represent this func- 
tion as 

G> (h,t 2 ) = -i (c,(*i)cj(t a )) = - l - {n |e^ /T c Q (ti)4(t 2 ) 

n 

= -^Y^e- E ^ T (n\c a \m)(m\cl\n}e^ E "- E ^- t2 K (246) 

nm 

In Fourie representation 

G>p(e) = -^^ e - £ " /T (n|c Q |m)(m|4|n) ( 5(£; n - E m + e). (247) 

nm 

Similarly, for the lesser function we find 

G<p(t) = ^^e-^/ T (n|4|m)(m|c Q |n) ( 5(ii;„ l -i?„ + e). (248) 

nm 

Now we can use these expressions to obtain some general properties of 
Green functions without explicit calculation of the matrix elements. Exchang- 
ing indices n and m in the expression (|248p and taking into account that 
E m = E n — e because of delta-function, we see that 

GUe) = -e-*' T G<p{e). (249) 

From this expression and relation (|236[) , which can be written as 

A af3 (e)=i[G> (e)-G< (e)] (250) 

we derive ()245l) . 



3.2.3 Interaction representation 

In the previous lectures we found that nonequilibrium Green functions can 
be quite easy calculated for free particles, and equations of motion for one- 
particle Green functions (the functions which are the averages of two creation- 
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annihilation operators) can be formulated if we add interactions and time- 
dependent perturbations, but these equations include high-order Green func- 
tions (the averages of three, four, and larger number of operators). The 
equations can be truncated and formulated in terms of one-particle Green 
functions in some simple approximations. However, systematic approach is 
needed to proceed with perturbation expansion and self-consistent methods 
(all together is known as diagrammatic approach) . The main idea of the dia- 
grammatic approach is to start from some " simple" Hamiltonian (usually for 
free particles) and, treating interactions and external fields as a perturbation, 
formulate perturbation expansion, and summarize all most important terms 
(diagrams) in all orders of perturbation theory. The result of such procedure 
gives, in principle, nonperturbative description (ordinary mean-field theory is 
the simplest example). The starting point of the method is so-called interac- 
tion representation. 

Let us consider the full Hamiltonian H as the sum of a free-particle time- 
independent part Hq and (possibly time-dependent) perturbation V(t) (note 
that this "perturbation" should not be necessarily small) 

H = H + V(t). (251) 

We define new operators in interaction representation by 

f'(t) = e iH tjS e -iH t^ ( 252 ) 

where f s is the time- independent Schrodinger operator. This is equivalent 
to the time-dependent Heisenberg operator, defined by the part Hq of the 
Hamiltonian. For a free-particle Hamiltonian Hq the operators f T {t) can be 
calculated exactly. 

A new wave function corresponding to (j252| is 

^(t) = e l " ot <P s {t). (253) 

It is easy to see that transformation (|252p . (|253p is unitary transformation 
and conserves the average value of any operator 

(W s \f 3 \W s ) = i&lf 1 ]* 1 ). (254) 

Substituting (]253[) into the ordinary Schrodinger equation, we derive the 
equation 

dW 1 - r J 
i— = V I (t)V I , (255) 

where V 1 (t) = e lH o t yS ^ e -iH t - 1S m ^ e intreraction representation. 

Equation ([255]) seems to be quite simple, however the operator nature of V 
makes this problem nontrivial. Indeed, consider a small time-step At. Then 
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<?(t + At) = [l - iV s (t)At\ #(t) = exp-^ S W At #(t), (256) 

linear in At term can be transformed into the exponent if we understand the 
exponential function of the operator in the usual way 



exp A = 1 + A + ^A 2 



(257) 



and assume that only linear term should be taken at At — > 0. 

If we now repeat this procedure at times U with step At, we obtain finally 



with 



t 

S(t,t ) = 11 exp (-iV 1 (ti)At 



(258) 
(259) 



this product, however, is not simply exp ^— i J V 1 (t')dt' j in the limit At — > 

0, because operators y(i') are not commutative at different times, and for 
two noncommutative operators A and B e A+B ^ e A e B . 

In the product (|259|) operators at earlier times should be applied first, 
before operators at later times. In the limit At — > we obtain 



(260) 



S(t,t Q ) = Texp ^-i j V'it^dt 

where T is the time-ordering operator ("-" for fermionic operators) 

ii(fi)B(ta) if h>t 2 , 



T 



(261) 



±B(t 2 )A(t 1 ) if ti < t 2 . 



Of cause, expression (|260p is defined only in the sense of expansion (|257D . 
Consider for example the second-order term in the time-ordered expansion. 



t -.2 

V^t^dt' 



= T 



V 7 (*')<**' / V'it^dt 



In 



dt' dt"V I (t')V 1 (t")+ dt" dt'V'it'^V 1 ^'). 

to J to J to J to 

(262) 



t" 



If we exchange t' and t" in the second integral, we see finally that 
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T 



t -|2 



V 1 (t')dt' 



to 



2 f dt' f df'^it'^it"). 

Jt Jt 



(263) 



(i) Properties of S(t,to) 

S is the unitary operator and 



S -1 (t,to) = (*,*(>) = Texp ^ijf t/ 7 



(264) 



where T is time-anti-ordering operator. Some other important properties are 

S-^Mo) = S(t ,t), (265) 
S(t3,t2)S(t 2 ,ti) = £(t3,ti), (266) 
S- 1 (t2,fi)S- 1 (i 3 ,f2) = S- 1 ^,^)- (267) 

Finally, we need the expression of a Heiscnbcrg operator, defined by the 
full Hamiltonian H = H + V(t), through an operator in the interaction rep- 
resentation. The transformation, corresponding to (|258|) . is given by 



f(t) 



_ „-iH to q-1 



S- 1 (f,i )/ J (t)S(f,to)e ifl °*°, 



(268) 



and the state ^ (to) is related to the Heisenberg time-independent wave 
function by 



I' 1 {t ) = e lHoto < T ' s 



(269) 



in accordance with our previous discussion of averaging we assume that at 
time t = to Heisenberg operators coincide with time-independent Schrodinger 
operators f H {to) = f s , and Schrodinger wave function coincides at the same 
time with Heisenberg time-independent wave function <P s (t ) = 'I fH . To avoid 
these additional exponents in (|268p we can redefine the transformation to the 
interaction representation as 



/*(*) 



iH (t-ta) fS e -iH (t-t ) 



(270) 



in accordance with the transformation (|206|) for time-independent Hamilto- 
nian. Previously we showed that free-particle Green functions are not de- 
pendent on to for equilibrium initial condition, if we want to consider some 
nontrivial initial conditions, it is easier to formulate these conditions directly 
for Green functions. Thus below we shall use relations 



f H (t) = S-'(t,t )f 1 (t)S(t,t ), 



(271) 



and 
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V 1 (t ) = <P S (t a ) = <P H . (272) 

(ii) Green functions in the interaction representation 

Consider, for example, the lesser function 

G< i8 (ti I ta)=*(4(t2)ca(ti)) = i(* H \c\(t 2 )c a (t 1 )\v H ), (273) 

c-operators here are Heisenberg operators and they should be replaced by 
operators c 1 (t) = c(t) in the interaction representation: 

G<,(tl,t 2 ) =i(^|5- 1 (t 2 ,to)4(*2)^2,to)5- 1 (tl,to)Ca(tl)5(t 1) t )|^). 

(274) 

Using properties of 5 operators, we rewrite this expression as 

G<p{t u t 2 ) = % (S(t , t 2 )c\(t 2 )S(t 2 , h^MS^to)) . (275) 

3.2.4 Schwinger-Keldysh time contour and contour functions 

(i) Closed time-path integration 

Now let us introduce one useful trick, so-called closed time-path contour of 
integration. First, note that the expression of the type 

f H (t) = 5- 1 (t,t )/ / (i)^(t,to) - Te 4/ 'o ^ (t ' )dt 7 7 (t)Te~ 4/ <o (276) 
can be written as 

f H (t) = T Ct exp V'^dt^ f J (t) 7 (277) 

where the integral is taken along closed time contour from to to t and then 
back from t to to 

/ dt'= dt' + / dt', (278) 
Jc t Jt Jt 

contour time-ordering operator Tc t works along the contour C-t, it means 
that for times t~ * it is usual time-ordering operator T, and for times f - it is 
anti-time-ordering operator T. Symbolically 

T Ct J dt' = T J dt'+f J dt'. (279) 
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Consider now the application of this closed time-path contour to calcu- 
lation of Green functions. It is convenient to start from the time-ordered 
function at t 2 > t\ 

(t (B(ta)A(ti)) ) = (S(t , t 2 )B{t 2 )S(t 2 ,t 1 )A(t 1 )S(t 1 ,to)) , (280) 

here A(t) and B(t) are Hciscnbcrg operators, A(t) and B(t) are operators in 
the interaction representation, in the case of fermionic operators the addi- 
tional minus should be added for any permutation of two operators. 

Using the properties of the S'-operator. we transform this expression as 

(5(*o,*a)-B(ta)^(t 2 ,*i)i(ti)5(ti,*o)) = (5- 1 (t 2 ,to)B(t2)5(t 2j ti)I(t 1 )5(t 1 ,t )) = 

= (s-\oc, t )S(ac, t 2 )B(t 2 )S(t 2 , ti)A(h)S(h, to)) = (s^T ^B(t 2 )A(h)s) ) , 

(281) 

where we defined operator 

S = S(oo,t ). (282) 
Using contour integration, it can be written as 

(T (B(t 2 )A(t x )) ) = (t c (§ c B(t?)A(t?j) ) , (283) 



S c = T c cxp (-% J V r (lf)dA , (284) 

contour C goes from to trough t\ and t 2 , and back to to - If t 2 > t\ it is obvious 
that contour ordering along C~* gives the terms from S(t\,to) to B(t 2 ) in 
(|280p . The integral over the back path gives 

T c exp (-i J V I (t/)dt'^j = Tcxp (-i ° V 1 ' (t')dt^j = 

= f cxp jf 2 V'^dt'^j = S-\t 2 ,t ) = S(t ,t 2 ). (285) 

For t 2 < ti the operators in (|280|) are reordered by T-operator and we 
again obtain (283]). 

The lesser and greater functions are not time-ordered and arguments of 
the operators are not affected by time-ordering operator. Nevertheless we 
can write such functions in the same form (|283p . The trick is to use one 
time argument from the forward contour and the other from the backward 
contour, for example 

(fl(ta)i(ti)) = (t c (s c B<$-)A(t?)) ) , (286) 
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here the time t\ is always before t 2 . 



(ii) Contour (contour-ordered) Green function 



Now we are able to define contour or contour- ordered Green function - the 
useful tool of Kcldysh diagrammatic technique. The definition is similar to 
the previous one 



G^(ri,r 2 ) = -i (t c (c Q (r n )4(T 2 )) ) , 



(287) 



where, however, n and r 2 are contour times. This function includes all 
nonequilibrium Green functions introduced before. Indeed, depending on con- 
tour position of times we obtain lesser, greater, or time-ordered functions 
(below we give different notations used in the literature) 



G a/3( T l,T- 2 ) = < 



Ti G G%r 2 G C~ 
Ti G C^,r 2 G G^ 
Ti,T 2 G : 



-i (Tc a (t 1 )c^(t2) 
-i /c a (ii)c^(i 2 ) 

-i (Tc Q (ti)c]j(t2) 



G— or G T {t u t 2 ), 
G+- or G>(t 1 ,t 2 ), 
G-+ or G<(t 1 ,t 2 ), 



=» G++ or G T {t u t 2 ). 
(288) 

These four functions are not independent, from definitions it follows that 



G< + G> 


= G T + G T , 


(289) 


and anti-hermitian relations 






Gap(tl,t 2 ) = 


-~G T p a (t 2 , ti), 


(290) 


Gap{tlM) = 


-G < p a {t 2 , ti), 


(291) 




p a {t 2 , t x ). 


(292) 



It is more convenient to use retarded and advanced functions instead of 
time-ordered functions. There is a number of ways to express G R and G A 
through above defined functions 



G R = 9(h - t 2 ) [G> - G<] =G T -G< 
G A = 9(t 2 - ti) [G< - G>] =G T -G> 



G> 
G< 



G T 1 
G f . 



(293) 
(294) 
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(iii) Contour Green function in the interaction representation 

In the interaction representation one should repeat the calculations performed 
before and given the expressions (|275j) . (|280|) . and then replace usual times 
by contour times t, so we obtain 

(295) 

Using contour integration, it can be written as 

G%{n,T 2 ) = -i (t c (ca(ri)ct ( T2 ))) = -i (t c (Sc*»(n)c£(7s))) , 

(296) 

S c = T c cxp (-i J V^t'jdt'^j . (297) 



3.3 Current through a nanosystem: 
Meir- Wingreen-Jauho formula 



Now we consider the central point of the NGF transport theory through 
nanosystcms - the Meir- Wingreen-Jauho current formula [5TJ [531 ISS] ■ This 
important expression shows that the current can be calculated, if the spectral 
and kinetic Green functions of the central system are known, and it is exact 
in the case of noninteracting leads. The details of the derivation can be found 
in the above cited papers, so we only briefly outline it. 



(i) Derivation by the NGF method 



In the absence of interactions in the leads (besides the tunneling) one can 
derive the following exact expression for the lead-system function: 



G< ika (e) =E^W [G%iz)9f k ^)+G< p {e)g? ka (e) 





(298) 



where ^^(e) and g^ a (e) are Green functions of isolated leads, Substituting 
it into (|198p . we obtain for the current 



2e 

Ji(t) = - 



de 

— Re 

2tt 



E *W^L,/3 [G% (e)9f k Ae) + G< (e)g? ka (e) 

k(7,a0 



(299) 
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For equilibrium right or left lead Green functions we obtain directly 

gUh -h)=i (cUhMh)) = ^f a A^)e- l{ ^ +eip){tl - t2 \ (300) 



= -iflfa -i 2 )e- i ( 6 *" +e ^( tl - ta ), 
(301) 

ff &(ti - t a ) = z0(t 2 - ti) ^ [c ka (h), cl(t 2 )] + ^ = i6(t 2 - iOe-^+^-H 

(302) 

or after the Fourier transform 

9kM = JdkAti - t 2 )e^ tl - t2) d(h - ta) = 2™/°(e^) ( 5(e - e fc(T - etp), 

(303) 

fl^(e) = -2tu[1 - f°(e ka )}6(e - e ka - etp), 



ffter(e) 



1 



e - efeo- - e<^ + zO ' 
1 

e - e ka - eip - iO ' 
/5(e) 



exp (^) + 1 



(304) 
(305) 

(306) 
(307) 



Using the level-width function (below without spin polarization of the 
leads) 



k(7 <T 

(308) 

and changing the momentum summation to the energy integration / p(e k )d 

k J 

we obtain the following expression for the current 

^Tr {r t (e - e Vi ) (G<(e) + /°(e - e Vi ) [G R (e) - G A (ejj ) } , 

(309) 

where /j is the equilibrium Fermi distribution function with chemical po- 
tential /Xj. Thus, we obtain the well-known Meir-Wingreen formula. Note, 
that we use explicitly the electrical potential of the leads in this expression. 
It is important to mention, that at finite voltage the arguments of the left 
and right level- width functions are changed in a different way, which means, 
in particular, that the known condition of proportional coupling F L = \r R 



Ji=L,R — 
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can be fulfilled only in the wide-band limit, when both functions are energy 
independent. 



(ii) Different forms of the MWJ formula 

In a stationary state J R — —Jl = J and one can use the symmetric form 



J= — 

2% I 2ir 



J ^r{[r L (e-e^ L )-r R (e-e m )]G<(e) 



+ [r L (e - e<p L )fi{e - e^ L ) - r R (e - ecp R )f° R (e - e<p R )] [c R (e) - G A (e)] } . 

(310) 

For the proportional coupling /i(e) = \r R (e) in linear response (ipi de- 
pendence of 2~j is ignored!) 

(311) 

A = i(G R —G A ) is the spectral function. This expression is valid for nonlinear 
response if the energy dependence of r can be neglected (wide band limit). 



(iii) Nonintcracting case 

Finally, in the nonintcracting case it is possible to obtain the usual Landaucr- 
Biittikicr formula with the transmission function 

T(e) = Tr [r L (e - e<p L )G R (e)r R (e - e m )G A (e)] . (312) 

This expression is equivalent to the one derived earlier by the single-particle 
Green function method. 

We should stress once more that this formula is valid for finite voltage. 
Therefore, the voltage dependence of the level-width functions is important. 



3-4 Nonequilibrium equation of motion method 

Now we start to consider the case of interacting nanosystems. Although the 
MWJ current formula is exact, the problem to find the Green functions of 
the central region is sometimes highly nontrivial. At the present time there 
are several techniques developed to solve this problem. 

Nonequilibrium equation of motion (NEOM) method is the simplest ap- 
proximate approach. In spite of its simplicity, it is very useful in many cases, 
and is very convenient for numerical implementation. In this section we con- 
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sider only a general formulation, some particular examples are considered 
further. 

We start from the general definition of a Green function as the average of 
two Heisenberg operators A(t) and B(t), denoted as 



A(ti),fl(i 2 ) 



R,A,< 



The particular definitions of the averages for spectral and kinetic functions 
are 



A(h),B(t 2 ) 



R 



-i6{h-t 2 )( A(h),B(t 2 ) 



(313) 



where upper sing here and below is for boson functions, lower sing for 
fcrmions, 

'AihlB^y = -i(i(tx),S(t 2 )). (314) 

The equations of motion for NGF arc obtained from the Heisenberg equa- 
tion of motion for operators 



dA 



A,H 



= AH - HA, 



(315) 



for any Heisenberg operator A{t). Here and below all Hamiltonians arc time- 
independent. We consider the stationary problem. 



(i) Spectral (retarded and advanced) functions 



Let us start from a retarded function 

R 



A(h),B(t 2 ) 



ie(t 1 -t 2 )( A(h),B(t 2 ] 



(316) 



Taking the time derivative we obtain 

.d_ 

''dh' 



^/(i(ii),B(i 2 )5) =S(t 1 -h)( A(ti),S(*i) 



A(Jti),H\_,B(t a ) 
(317) 



where the first term originates from the time-derivative of the (9-function, and 
the equation (|315p is used in the second term. 

In the stationary case the Fourier transform can be used 



(e + irj)({A,B 



A, B 



A, H 



,B 



(318) 



Now let us assume that the Hamiltonian can be divided into " free particle" 
and "interaction" parts H = Hq + Hi, and [A, Hq]^ = IqA. (The simple 
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example. For the free particle Hamiltonian Hq = epdLdp and the operator 

A = d) a one has [A,H ]- = Y,p e /3^L dpdp]- = e a d+, e = e a is simply a 
number. In general, eo is some time-independent operator). So that 



(e + ir)-e )({A,B 



A,B 



A, ^ 



B 



(319) 



the second term includes interaction and can not be easy simplified. 

It is convenient now to introduce the "free particle" function gf as a 
solution of the equation 



(e + irj - e )<?f = 1. 



(320) 



Now we multiply the right and left parts of (|319[) by gf . Using the function 
g R (t) = J g R e~ let ^ we can write the time-dependent solution of (|317[) as 



A(h),B(t 2 ))) -ta)< A(h),B(h) 



+ 9 R {U-t')(( 1(0,^1 ,B(t 2 )]) dt'. (321) 



(ii) EOM on the Schwingcr-Keldysh contour 

The calculation of the lesser functions by the EOM technique requires some 
care. To demonstrate it let us compare the EOM for retarded and lesser 
functions of free particles. 

The equation for g R p is (assuming the diagonal matrix e a p) 

(e + iv-e a )g^ = 5 a p, (322) 

from which the free-particle Green function is easily obtained. 
At the same time for the lesser function we have the equation 

5^ = 0, (323) 

from which, however, the free-particle lesser function = 27r/o(e)<5(e — 
£a)fia(3 can not be obtained. 

The problem can be generally resolved by using the EOM on the Schwinger- 
Keldysh time contour. Contour-ordered Green function is defined as 

((i(n),.S(7*))) C = -*(T c (i(r 1 ),B(r 2 ))), (324) 

where A(n) and Bfa) are two Heisenberg operators, defined along the con- 
tour. 

Taking the time derivative we obtain the equation 
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A( Ti ),H^_,B(t 2 ) 
(325) 



in the stationary case this equation can be formally solved if one applies 
the Fourier transform along the contour, or perturbation expansion in the 
interaction representation (Niu et al. 1999). Using the free particle solution 
g c (ti — t 2 ) we can write the time-dependent solution as 



A( Ti ),B(t 2 ) 



r(n-r 2 )( AinlBin) 



(ri-r')(( A(/),&! ,B(r 2 ))) dr'. (326) 



(iii) Kinetic (lesser) function 

Applying now the Langreth rules (see the next section for details), which 
shows, that from 



C(r 1 ,r 2 )= / A(T 1 ,T 3 )B{T 3 ,T 2 )dT 3 
Jc 



(327) 



it follows 



C R (h,t 2 ) = f A R {t u t 3 )B R {tz,t 2 )dh, (328) 
C<(h,t 2 )= J (A R {t 1 ,t 3 )B R {t 3: t 2 )+A<(t ll t 3 )B A {t 3 ,t 2 ))dt 3 , (329) 

we get ()321j) for the retarded function, and 



i4(ti),s(ta)» =g < (h-t 2 )( 4(ti),s(ti) 







+ / 9 R (h - 










+ J ^(h- 




i(i')X 


for the lesser function. And the Fourier transform is 




A, B 


>»•"« 


A, Hi 


]..*>:■ 



A,H t 



,B 



(331) 
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3.5 Kadanoff-Baym-Keldysh method 
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Now we review briefly the other approach. Kadanoff-Baym-Keldysh (KBK) 
method systematically extends the equilibrium many-body theory to the 
nonequilibrium case. Potentially, it is the most powerful approach. Below 
we give a simple introduction into the method, which is currently actively 
developed. 



(i) Perturbation expansion and diagrammatic rules for contour functions 

We found that Green functions can be written in the interaction represen- 
tation with a help of the S'-operator. For example, time-ordered fermionic 
Green function is 

G^/j(*i > t2)=-*(r(c a (ti)cJ(ta))) = -i(5- 1 T(c a (t 1 )c^(t 2 )5)), (332) 
using "usual" S-operator 

S = S(oo,t ) =Texp y-ij^ V^t^dt'j , (333) 

or 

Glp(h,t 2 ) = -i (Tc (ca(t?)cl(t?)§c) ) , (334) 
using "contour" Sc-operator 

S c = T c exp f-i J V T (t')dA . (335) 

We first consider the zero temperature case, when one can set to = — oo, 

S = 5(oo,-oo) =Tcxp (-i J V^t^dt^j, (336) 

and assume that interaction is switched on and switched off at t — > +oo 
adiabatically. This condition is necessary to prevent excitation of the system 
from its ground state. The other necessary condition is that the perturbation 
is time-independent in the Schrodingcr representation. In this case if the 
initial state \\P(t = — oo)) = I'f'o) is the ground state (of free particles), then 
the final state \&{t = +oo)) = S\& a ) = e i9 \&°) is also the ground state, only 
the phase can be changed. Now, using the average value of the S'-operator 

(S) = (V°\S\V°) = e w (V°\V°) = e l \ (337) 



we obtain 
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S\V°) = (S)\&°), (338) 

and 



(^"l^T 1 = ^J-. (339) 
V ' (5) 



So that (|332[) can be written as 

(c a {ti)cl(t 2 )S 
~W) 



T 

G^/j( f i>*2) = - i — — — ■ ( 340 ) 



Now we can expand the exponent (note that S'-operator is defined only in 
the sense of this expansion) 

S =Tcxp -i / V\t')dt> = TJ2 { —r / <- / < 

V J-oo / n=0 n - J-oo J-oo 

(341) 

and numerator and denominator of the expression (|340|) arc 

(T(Uh)cl(t 2 )s)) = f2^ r dt[... r dt'/Tc a ( tl )di(t 2 )r(t[)...r(t' n )), 

' ' _ 71. J — oo J —oo x ' 

(342) 



n=0 



( S ) = Y, -J-J J^-j x < (TV I (t[)...V I (0) . (343) 

These expressions are used to produce the perturbation series. 
The main quantity to be calculated is the contour Green function 

G(l,2) = G%(n,r 2 ) = -i(T c (ca(n)cj(7*))) , (344) 

where T\ and t 2 are contour times. Here l c = a,T\. 

The general diagrammatic rules for contour Green functions are exactly 
the same as in the usual zero-temperature technique (we call it standard 
rules). The correspondence between diagrams and analytical expressions is 
established in the following way. 

1. Open bare electron line is iGo(l, 2). 

2. Closed bare electron line is no(l) = n>a (ti). 

3. Bare interaction line is — iv(l, 2). 

4. Self-energy is —iS{l, 2). 

5. Integration over internal vertices, and other standard rules. 



Green function techniques in the treatment of quantum transport 73 

(ii) Langreth rules 

Although the basic equations and diagrammatic rules are formulated for con- 
tour Green functions, the solution of these equation and final results are much 
more transparent when represented by real-time spectral and kinetic func- 
tions. 

As in the ordinary diagrammatic technique, the important role is played 
by the integration (summation) over space and contour-time arguments of 
Green functions, which is denoted as 

[dl c = J2 [ dn. (345) 
J a Jc 

After application of the Langreth rules, for real-time functions these integrals 
become 

/POO 
dl=^2 dtl - ( 346 ) 

The Langreth rules show, for example, that from 

C(n,r 2 )= / A(r u r 3 )B(T 3 ,T 2 )dT 3 (347) 
Jc 

it follows 

C R (h, t 2 ) = J A R (t 1 ,t 3 )B R (t 3 , t 2 )dt 3 , (348) 
C<(h,t 2 )= J {A R (t ll t 3 )B<(t 3 ,t 2 ) + A<(t 1 ,t 3 )B A (t 3 ,t 2 ))dt 3 . (349) 

The other important rules are: from 

C{t 1 ,t 2 ) = A(t 1 ,t 2 )B(t 1 ,t 2 ) (350) 

it follows 

C R : (t ! , t 2 ) = A R (t x , t 2 ) B R \t x , t 2 ) + A R (h , t a )B< (t 1 , t 2 ) + A< (t 1 , t 2 )B R \t x (mi) 
C<(h,t 2 ) = A<(t 1 ,t 2 )B<(t 1 ,t 2 ), (352) 

and from 

a(ri,r 2 ) = A(r 1 ,T 2 )5(r 2 ,Tx) (353) 

it follows 

C R ih,t 2 ) = A R (t 1 ,t 2 )B<(t 2 ,t 1 ) + A<(t 1 ,t 2 )B A (t 2 ,t 1 ), (354) 
C<{t u t 2 ) = A<{tx,t 2 )B>(t 2 ,tx). (355) 
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Fig. 15 Diagrammatic representation of the first-order self-energy. 



(iii) First-order sclf-cncrgy and polarization operator 

Consider, as an example, the first order expression for the self-energy, shown 
in Fig. 1151 Following the diagrammatic rules, we find 

Si(l, 2) = 5{1 -2) J v(l, 3)n (3)d3 + iv(l, 2)G (1, 2), (356) 

where the first term is the Hartree contribution, which can be included into 
the unperturbed Green function Go (1,2). This expression is actually sym- 
bolic, and translation from contour (Keldysh-time) to real-time functions is 
necessary. Using the Langrcth rules, one obtains 

E?(l,2) =5(1+ -2) J w fl (l > 3)n (3,3)d3 + it; H (l > 2)G?(l > 2) 

+ iv<(l, 2)G*(1, 2) + iv R (l, 2)G<(1, 2), (357) 

E<(l,2) = iv<(l,2)G<(l,2). (358) 

There is no Hartree term for lesser function, because the times T\ and t-i are 
always at the different branches of the Keldysh contour, and the <5-function 
S(ti — t-i) is zero. 

In the stationary case and using explicit matrix indices, we have, finally 
(r = t\ — 1<£.) not to mix with the Keldysh time) 

^3 1) (r)=^(r + )^E 7 <(0)4 0) 
+™%{t)G I $ ) t) +iv< p {r)Glf ) (r) +rv^{r)G < a f\r) 1 (359) 

^ 1 \r)=iv< (r)G< O) (r), (360) 
and we define the Fourier transform of the bare interaction 

S2r(0) - / v*(T)dT. (361) 



Finally, the Fourier transforms are 
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Fig. 16 Diagrammatic representation of the first-order polarization operator. 




(0) 



2^ 



(362) 



(363) 



The second important function is the polarization operator ("self-energy 
for interaction"), showing in Fig. 1161 Following the diagrammatic rules, we 
find 

211(1,2) = -iGo(l,2)Go(2,l), 



note the order of times in this expression. 
Using the Langrcth rules, 

n* (1, 2) = iG$(l, 2)G<(2, 1) + 1 G<(1, 2)G#(2, 1) 

JT<(l,2) = iG<(l,2)G>(2,l). 
And in the stationary case, restoring the matrix indices 



n 



at, ( £ ) = -* Gj(r)G^(-r)+G a 7(r)G^(-r) 



Ag 1) W = -i^(r)G>i«(-r). 



(364) 

(365) 
(366) 

(367) 
(368) 



In the Fourier representation 



n 



Ctj3 



(r) = -* / £ G^\e')G<°\e> - e) + G<J\e')G^\e' - e) 



2vr 



(369) 

n$Hr) =-iJ ^G<°\e')G>i°\e' - e). (370) 

These expressions are quite general and can be used for both electron- 
electron and electron-vibron interaction. 

For Coulomb interaction the bare interaction is is v(l, 2) = U a p8{Ti — T2), 
so that 
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Fig. 17 Diagrammatic representation of the Dyson equations. 



v R (l,2) = U a pS(t+-t 2 ), (371) 
u<(l,2) = 0. (372) 



(iv) Self-consistent equations 

The diagrams can be partially summed in all orders of perturbation theory. 
The resulting equations are known as Dyson equations for the dressed Green 
function G(l,2) and the effective interaction W(l,2) (Fig. [T7]) . Analytically 
these equations are written as 

G(l, 2) = G (l, 2) + / / G (l, 3)17(3, 4)G(4, 2)d3d4, (373) 



W(l ) 2)=u(l,2)+ // v(l, 3)J7(3,4)W(4,2)d3d4. (374) 



In the perturbative approach the first order (or higher order) expressions 
for the self-energy and the polarization operator are used. The other pos- 
sibility is to summarize further the diagrams and obtain the self-consistent 
approximations fFigs. fT8H9p . which include, however, a new unknown func- 




Fig. 18 Diagrammatic representation of the full self-energy. 



®- 




Fig. 19 Diagrammatic representation of the full polarization operator. 
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• = • + 




Fig. 20 Diagrammatic representation of the vertex function. 



tion, called vertex function. We shall write these expressions analytically, 
including the Hartree-Fock part into unperturbed Green function Gq(1,2). 



The equation for the vertex function can not be closed diagrammatically 
(Fig. [20]) . Nevertheless, it is possible to write close set of equations (Hedin's 
equations) , which are exact equations for full Green functions written through 
a functional derivative. Hedin's equations are equations (|373p - (|376[) and the 
equation for the vertex function 



r(l;2,3) = <J(l,2)(J(l,3)+ / / / / G(4, 6)G(7, 5)r(l; 6, 7) ' ( d4d5d6d7. 




(375) 




(376) 




(377) 
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4 Applications 

4-1 Coulomb blockade 

In Section 2 we have seen that Coulomb blockade phenomena mediated 
by electron-electron interactions on a quantum dot can be dealt with in 
a straightforward way by using master equation (ME) approaches, which 
are based on Fermi's Golden Rule. PUJ UM OH HH E22 EH UTS 
However, due to its intrinsic perturbative character in the lead-dot coupling, 
ME techniques cannot cover the whole interaction range from weak-coupling 
(Coulomb blockade), intermediate coupling (Kondo physics), up to strong 
coupling (Fabry- Perot physics). It is thus of methodological and practical 
interest to develop schemes which allow, in a systematic way, to describe the 
three mentioned regimes also in out-of-cquilibrium situations. As stated in 
the introduction, we believe that Green function techniques are such a tool; 
in this section we will show how a non-equilibrium treatment of the Hubbard- 
Anderson model together with appropriate approximations allow us to repro- 
duce the well-known Coulomb blockade stability diagrams obtained with the 
master equation approach (see also Section 2). For the sake of simplicity we 
will deal with the problem of single and double-site dots in the CB regime, al- 
though the method can be straightforwardly extended to multi-level systems. 
Our purpose is to study the problem of a two site donor/acceptor molecule in 
the CB regime within the NEGF as a first step to deal with the phenomenol- 
ogy of a rigid multilevel island. The nuclear dynamics (vibrations) always 
present in molecular junctions could be then modularly included in this the- 
ory. Our method can be calibrated on the well-studied double quantum dot 
problem |176[ 1192] and could be possibly integrated in the density functional 
theory based approaches to molecular conductance. The Kondo regime would 
require a separate treatment involving more complex decoupling schemes and 
will be thus left out of this review, for some new results see Ref. [213| (EOM 
method) and Refs. 224 ( 225, 1226j (the self-consistent GW approximation). 

The linear conductance properties of a single site junction (SSJ) with 
Coulomb interactions (Anderson impurity model), have been extensively 
studied by means of the EOM approach in the cases related to CB [203 , 204 
and the Kondo effect. |205j Later the same method was applied to some 
two-site models. [2061 [2071 [2081 [214] Multi-level systems were started to be 
considered only recently. [210[ 1211] For out-of-equilibrium situations (finite 
applied bias), there are some methodological unclarified issues for calculating 
correlation functions using EOM techniques. [212[ I213[ 1214] We have devel- 
oped an EOM-based method which allows to deal with the finite-bias case in 
a self-consistent way. [209 
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4.1.1 Nonequilibrium EOM formalism 
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(i) The Anderson-Hubbard Hamiltonian 



We consider the following model Hamiltonian (which can be called the multi- 
level Anderson impurity model, the Hubbard model, or the quantum cluster 
model) 

H = e Q!| a<4d l a+^ E ^a/^a^/a+E hkac\ ka c iklJ + E (v. lklJ , a c\ kcr d a + h.c^j 

a0 a0 ika ika, a 

(378) 

electrical potentials are included into the energies <Li ka = ti ka + etpi(t) and 

This model is quite universal, describing a variety of correlated electron 
systems coupled to the leads: the Anderson impurity model, the multilevel 
quantum dot with diagonal nonintcracting Hamiltonian e a p, a system (clus- 
ter) of several quantum dots, when the off-diagonal matrix elements of e a( g 
describe hopping between individual dots, and, finally, the ID and 2D quan- 
tum point contacts. 



(ii) EOM for Heisenbcrg operators 

Using the Hamiltonian (|378|) one derives 
. dc ika 



dt 



Cikcr , H 



^ika^iko 



E^cA, (379) 



^ = -hkAka ~ E t&r,«4» (380) 



i ~Q~f = E l °-f id l } + E U ^pnpd a + E V <Lk*, a c ika, (381) 

0=£a ika 

^ = - E *«f>4 - E - E w4,> (382) 

0=£a ika 



ika 

+ E £~f0 d *j d i3 ~ E ^7 rf l rf 7- (383) 

a 

These equations look like a set of ordinary differential equations, but are, 
in fact, much more complex. The first reason is, that there are the equations 
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for operators, and special algebra should be used to solve it. Secondly, the 
number of Cjfc CT operators is infinite! Because of that, the above equations 
are not all sufficient, but are widely used to obtain the equations for Green 
functions. 



(iii) Spectral (retarded and advanced) functions 

Now we follow the general NEOM method described in the Section 3. Using 
(|38ip . we get the equation for G^ g = — i / \d a , 4 J 

(e + %ri)G*p - ]T ~^G^ = 5 a0 + J2 U ay G®% + ]T V^G*^ (384) 

7 7^a ika 

which includes two new functions: G^l R g and G^ ka p. 

The equation for Gf ka g is closed (includes only the function G R S intro- 
duced before) 

(e + hi - hkc)Gf k<Jt g = Vik^Gfg. (385) 

8 

The equation for 

G£$(*i-*a) = -W(ti-*a)^[d (*i)n 7 (*i),4(*a)] + ) 

is more complicated 

(e + iv)G^ R p -^2e aS Gf^ = n 7 S a0 + (6 aa - p a g)6 ai 
s 



V*k<r, 7 ^d a d\cika; 4}) ~Y Vik °<-y (( d « c lfe<7 d 7; 4 



+^ M ((d a d\ds; 4 )) - ^ Q 7 (( <4 4 ^7 i 4 )) ■ ( 386 ) 



R 



The equation (|386[) is not closed again and produces new Green functions 
of higher order. And so on. These sequence of equations can not be closed in 
the general case and should be truncated at some point. Below we consider 
some possible approximations. The other important point is, that average 
populations and lesser Green functions should be calculated self-consistently. 
In equilibrium (linear response) these functions are easy related to the spec- 
tral functions. But at finite voltage it should be calculated independently. 
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(iv) Kinetic (lesser) function 

Following the same way, as for the retarded functions (using only the defi- 
nitions of NGF and Heisenberg equations of motion) one derives instead of 
pg]) - p5]) 

e< ^Q/3 ~ X/ ^ a 1^lP = X ^ a ^a-u0 + X] Vika,aGfka,i3' (387) 
7 77^^ 

(e - hka)G<^ = £ Vi^sGfp, (388) 

+ X! ^L.a ((cifeo-%; 4)) + H ^^,7 ((^4^°-; 4}) 

- X! VW, 7 (^d a c\ ka d 7 ; 4}} 
+ X e 76 (^dad^dg; 4}) < -l]e57 ((^K; 4)^ ■ ( 389 ) 

To find Gf kl7 p we should divide the right parts by (e — e^), which is 
not well defined at e = e^. In the section 3 we considered the general 
prescription to avoid this problem, we use the equation (|331[) , and instead of 
(|555|) wc obtain 

Gfk*,f3 = 9L E V ^sG< 3 + 9 < a E Vtt ffl4 G^. (390) 
s s 

The equations (|387[) and (|389[) can be used without modifications because 
they include the imaginary parts (dissipation) from the lead terms. 

At this point we stop the general consideration, and introduce a powerful 
Ansatz for the NGF which is related both to the equation-of-motion (EOM) 
method and to the Dyson equation approach. [209] From the knowledge of 
the Green functions wc then calculate the transport observables. For clarity, 
we first describe our method in the more familiar problem of a single site 
junction, which is the well-known Anderson impurity model. Then wc apply 
it to a double quantum dot. The equations obtained below by the heuristic 
mapping method can be obtained straightforward from the general NEOM 
equations derived in this section using the same approximations as in the 
mapping method. 
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4.1.2 Anderson impurity model (single site) 



The Anderson impurity model is used to describe the Coulomb interaction 
on a single site: 

H = Hn + y](H a + H a v), 

where 

a 
k.a 

^aD = {^a,k,a C a,k,a d a + Va,k,v d a C a,k,a) ' 

k.a 

where d and c are the operators for electrons on the dot and on the left 
(a = L) and the right (a = R) lead, U is the Coulomb interaction parameter, 
e a is the a level of the quantum dot, while e£ a is the spin a level of lead a in k 
space, a =t, j. With the help of the EOM and the truncation approximation, 
we can get a closed set of equations for the retarded and advanced GFs 

(w - e„ - K /a )Gl{r = <W + UG^'\ (391a) 
(cu-e„-U- E*J*)GM T >* = {n 9 )S afT , (391b) 

where G Z J* = {{d a \d\)Y^, G^ r/a = «n^ CT |4» r / a and 

= + = £ ~ z¥^+ (392) 
are the electron self-energies. 

(i) Mapping on retarded Green functions 

For retarded GFs, from the EOM method, and with the help of Eqs. (|391a|l 
and (|391b]) . we can get 

G r = G T + G T UG^ 2)l = G + Go^ BOM G (1)r , 

where G r is single-particle GF matrix 
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fG\ T G r T A 

and G^r = Gi- 2 r / ' (n s ) . G r describes the single-particle spectrum with- 
out Coulomb interaction, but including the effects from the electrodes. 
£^® M = U {n s ) is the Hartree-like self-energy of our model. Since there 
is only Coulomb interaction on the site with the levels e a , the Fock-like self- 
energy is vanishing. 

Alternatively, by means of the Dyson equation and the second-order trun- 
cation approximation, taking Hartree-like self-energies = U(n s ) (= 
£p° M ), we can also get the retarded GFs as follows 

G r = Gl + GlS n G\, (393) 

where G\ = GJ, + Gq17 h Gq is the first-order truncation GF. 

Within the level of the second-order truncation approximation, we see that 
there is a map between the EOM results and the Dyson results: 

G r = G + G T S H G (1)r (EOM), (394a) 
I t 

G r = G + G T S u G\ (Dyson). (394b) 

Eqs. (|394|) prompts a way to include further many-particle effects into the 
Dyson equation, Eq. (|394b|) . by replacing the Dyson-first-order retarded 
Green function G\ with the EOM G^ r . Then one obtains already the correct 
results to describe CB while keeping the Hartree-like self-energy. 

(ii) Mapping on contour and lesser Green functions 

Introducing now the contour GF G, we can get the Dyson equation as fol- 
lows [73 El El E5] 

G = G + G SG, (395) 

where £ is the self-energy matrix. [66] 

According to the approximation for the retarded GF in Eq. (|393j) , we take 
the second-order truncation on Eq. ()395|) . and then get 

G = Go ~\~ GqX/^Gi : 

where Gi = Go + Go^ H Go is the first-order contour GF, and Go has already 
included the lead broadening effects. 

Similar to the mapping in Eq. (|394|) . we perform an Ansatz consisting in 

substituting the Dyson-first-order G[ /a/< with the EOM one G (1)r/a/< to 
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consider more many-particle correlations, while the EOM self-energy is used 
for the Dyson equation for consistency: 

G = Go + Go £ H Gi (Dyson), 

I T (396) 

G (EOM). 

Then, using the Langreth theorem [66] we get the lesser GF, 

G< = Gq + G Z ,H ^'G (1)< + G£2? H ' a G* 1 > a 

= G< + G C/G (2)< + G<J7G (2)a (397) 

where G^^ are GFs for U = 0, &u£ including the lead broadening effects, 



i.e. 



G< = g< + glE<Gl + g<Z a G% + g^G< 
Go /a = .9o /a + .9o /a ^ r/a G /a , 



with a/ a//< the free electron GF, and 



/ yWa/< 

r r/a/< = 



^I /a/ 7 ' 

^=iE a ^a/aH, and r a = i{Z T a -E*), f a (w) = f(w-fi a ), f is the cqui- 
librium Fermi function and \x a is the electro-chemical potential in lead a; 
Eu & are the retarded/advanced electron self-energies from Eq. (|392[) and 
G^/t = Gajr / { n s)> Performing the same Ansatz on the double- 
particle GF, from Eq. ([391b)) we can get 

G (2)< = G ,(2)r i: (2)< G (2)a 5 (393) 

with Si 2)< =S</{n 9 ). 

The lesser GFs in Eq. (|397j) can also be obtained directly from the general 
formula |66| 

— L^o 1 ^0 5 

with the help of the Ansatz in Eq. ((396)) . It should be noted that Eq. ([397) 
is very different from the lesser GF formula, 

G< = G r £<G a , (399) 

with the self-energy S < containing only contributions from the electrodes. 
The equation (|399[) is widely used for both first-principle 255, 236, 258J and 
model Hamiltonian calculations. |207j 
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Fig. 21 (Color) The stability diagram of a SSJ with e a = 2.0 eV, U = 4.0 eV, 
Fi, = = 0.05 eV. (a) The uncorrect result obtained by means of the widely used for- 
mula in Eq. (1 399 P for the lesser GF is not symmetric for levels e a and e CT + U. (b) Results 
obtained by means of our Ansatz in Eq. l|397( shows correctly symmetric for levels e CT and 



The numerical calculation results of conductance dependence on the bias 
and gate voltages by the two different NGF Eqs. (|397p and (|399|) are shown 
in Fig. [2D As we can see in the left panel, the adoption of Eq. p99[) results in 
an incorrectly symmetry-breaking in the gate potential. This wrong behavior 
is corrected in the right panel where Eq. (|397[) has been used. 

Note, that the expressions for the retarded and lesser functions, described 
above, can be obtained in a more formal way by the EOM method formulated 
on the Keldysh contour. 



(iii) Comparison with the master equation result 

In the single site model with two (spin-up and spin-down) levels it is possible 
to make the direct comparison between our Ansatz and the master equa- 
tion methods. For the latter, we used the well known master equations for 
quantum dots [TglHTgD] . 

In the Fig. [22] the typical curves of the differential conductance as a func- 
tion of the bias voltage at fixed gate voltage obtained by the two methods 
are shown together: there is basically no difference in the results obtained by 
these two methods. In the Fig. [23] the contour plot of the differential conduc- 
tance obtained by our Ansatz is shown. We do not present here the contour 
plot obtained by the master equation method because it looks exactly the 
same. 

It is quite clear from the presented figures that our Ansatz and the master 
equation method give essentially the same results in the limit of weak coupling 
to the leads. The systematic investigation of the deviations between the two 
methods at stronger tunneling will be presented in a separate publication. 
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Fig. 22 (Color) The comparison of the master equation method and our Ansatz for 
the differential conductance of the two level model with e-t = —0.35 eV, ei = —0.65 eV, 
U = 1.0 eV, V g = 1.0 V, T L = T R = 0.05 eV. 



It is important that our Ansatz can be applied straightforwardly to multi- 
level systems in the case when the exact eigenstates of an isolated system are 
unknown and the usage of the master equation method is not easy. In this 
paper we consider the simplest example of such a system, namely a double 
site case. 



&{2e 2 /h) 




r 00 



Fig. 23 (Color) The stability diagram (the contour plot of the differential conductance) 
calculated by our Ansatz for the two level model with parameters as in Fig. l22l The latter 
is indicated with a dash line at V K = 1.0 V. 
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4.1.3 Double quantum dot (two sites) 

We now return to the investigation of the DS J system (Fig. |24|) with Coulomb 
interaction on each site. The Hamiltonian is expressed as follows, 



where 



H = H D + H t + J2( H » + H <*v), 



i,(T 

H =Ve (a) c t c 

k,a 

H a D,a = ya,k,a C a,k,a d i,cr + ^a,k,cr d i,a C a,k 



k.a 



with i,j = 1,2 indicate the site, t is the constant for electron hopping between 
different sites. 

With the help of the EOM, and by means of the truncation approximation 
on the double-particle GFs, we obtain the closed form for the retarded GFs 
as follows 

(u - V - Sl a )G^ T = 5.^ + UtG?^ + t Gfc$, (400a) 
(<" - % a -Ui- 2^G<?W r = (n^S^ + t n if9 G%%, (400b) 
where the DSJ retarded GFs arc defined as 




Fig. 24 (Color) The general configuration of a double site junction. The levels ei,2 with 
charging energies U\ 2 are connected via t and coupled to the electrodes via the linewidth 
injection rates 7^. 
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G 



(2)([/,t)r 



(401) 

(402) 



Here i means 'NOT i\ and S\ a are the electron self-energy from leads. 

From Eqs. (|400ap . (|400b[l and performing the same Ansatz as in the case 
of SSJ, we can obtain the DSJ lesser GFs with Coulomb-interaction effects 
as follows 

G {u ' t]< (Lu) = (1 + G {u ^Sl)G {u)< {\ + £?G (u ' t)a ) + G iu ' t}r S<G (u ' t} lm) 
with 



(0 t 0\ 
i 
t 

\0 0t0j 

























°\ 


u = 











ih 






r — 






7* 





ik 
















c/ 2 y 












7*/ 



and 27 t < = 0. G^< is the DSJ lesser GF with the same form as Eq. (f397|) . 
but taking 



(404) 



where 7* indicates the line width function of lead a to site i, and Ui is the 
charging energy at site i. G r / a and G' 2 - )r / a are the GF matrix from Eqs. (|400a|l 
and (|400b[l . Here, in order to distinguish different GFs, we introduce the 
subscript '(U,t)' for the one with both Coulomb interaction U and inter-site 
hopping t, while l (U)' for the one only with Coulomb interaction. 

For our models, the lesser GFs in Eq. ([557]) . ([555)) and (|4"05]) . which are 
obtained with help of our Ansatz, can also be obtained by the EOM NEGF 
formula (|33 1 1) within the same truncation approximation. 

The current can be generally written as |81j 



r R )G^< + [/i»r L - h(Lo)r R ](G^ - g^)}, 



where the lesser GF is given by Eq. (I403p . The differential conductance is 
defined as 

dJ 



dVbias ' 

where the bias voltage is defined as Vbias = (Mr — MiJ/e. 
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Fig. 25 (Color) The stability diagram of a serial DSJ with ei i<T = £2,0- = —0.15 eV, 
Ui = Ui = 0.3 eV, t = 0.05 eV, 71=72= 0.02 cV, 7 £ = 7^ = ,V bias = 0.005V. The 
maximums of conductance are observed when the levels of the first site (ei, CT or ei jtT + U) 
are overlapped with the levels of the second site (£2,0- or e2,<r + U), and with the Fermi 
energy in the leads. The splitting of the four maximums is due to the hopping between the 
dots. 

(i) Serial configuration 

By taking l\ — = 0, we obtain a serial DSJ, which could describe 
the kind of molecular quantum junctions like the ones studied in Ref. [IB] . 
First, at small bias voltages, the conductance with the two gate voltages V g 
and V s was calculated, and the relative stability diagram was obtained as 
shown in Fig. 1251 Because of the double degeneracy (spin-up and spin-down) 
considered for each site and electrons hopping between the dots, there are 
eight resonance-tunnelling regions. This result is consistent with the master- 
equation approach. |176| 

Further, we studied the noncquilibrium current for large bias-voltages 
(Fig. [26]) . Because ei j<r and e2,cr are taken as asymmetric, for the case with- 
out Coulomb interaction, the I-V curve is asymmetric for ±Vbias, and there 
are one step and one maximum for the current. The step contributes to one 
peak for the conductance. When we introduce the Coulomb interaction to 
the system, the one conductance peak is split into several: two peaks, one 
pseudo-peak and one dip, while the current maximum comes to be double 
split (sec Fig. [26l) . The origin of this is in the effective splitting of the degen- 
erate level when one of the spin states is occupied and the other is empty. 
When both spin states are occupied, the degeneracy is restored. 

This process can be illustrated by the help of Fig. 1271 At zero bias- voltage, 
£2.0- is occupied and ei jCr is empty. Then we start to increase the bias voltage, 
a) The level £.2,0 + U is first opened for transport. It will contribute the 
first peak for conductance, b) Further, the levels e 2 ,er and e l cr come into the 
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Fig. 26 (Color) Current and conductance vs. bias-voltage of a DSJ far from equilibrium 
with parameters ei, a = 0.5 eV, e 2 , CT = -0.5 eV, Ui = U 2 = U = 0.2 eV, t = 0.07 cV, 
7l = Tr = 003 eV > Vg 2 = -Vgi = V hia a/4 and V K = -V h = V hias /2. The red curve 
represents the current, while the blue the conductance. The inset is the blow-up for the 
conductance peak split. The dash and dot-dash curves are for current and conductance 
with (7 = 0, respectively. 



transport window between the left and the right Fermi levels, resulting in 
the second peak, c) When the level + U comes into play, only a pseudo- 
peak appears. This is because there is only a little possibility for electrons 
to occupy the level e\ tCr under positive bias voltage, d) Levels £2,0- + U and 
ei.o- meet, which results in electron resonant-tunnelling and leads to the first 
maximum of the current. Then a new level ei i(T + U appears over the occupied 
ei.o- due to the Coulomb interaction, e) The meeting of 62.0- and ei j(T results 
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Fig. 27 (Color) The processes involved in the transport characteristics in figure l26l ei = 
ei jCT , e2 = £2,(ti The red line indicates electron resonant-tunnelling, a) The first conductance 
peak, b) The second conductance peak, c) The pseudo-peak of conductance, d) The first 
current maximum, and the red line indicates resonant tunnelling of electrons, e) The second 
current maximum for electron resonant tunnelling, f) The dip of conductance. 
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in electron resonant tunnelling. It means that t\_ a will be occupied, which 
leads to the appearance of a new level + U. Then £2,0- + U meets e\ M + U 
and another resonant tunnelling channel is opened for electrons. The two 
channels result in the second current maximum, f) finally, the level ei j0 - + U 
disappears if the level e\ M is empty. This means that a dip appears in the 
conductance. 

It should be noted that the characteristics of serial DSJ in Fig. [26] have 
showed some reasonable similarities to experiments of a single-molecule 
diode. [18] 

(ii) Parallel configuration 

If on the other hand, the two sites are symmetrically connected to the elec- 
trodes, possibly with a small inter-dot hopping, but with charging energies 
XJ\ and U2 fixed to different scales for transport. The resulting stability dia- 
gram contains both interference effects for parallel pathways and an overlap 
of U± and U2 stability diagrams, which we refer to a nesting characteristic, 
(see Fig. US). 

The physics of the weak lines in the figure can be understood by the 
help of charging effects. For simplicity, here we would ignore the site index 
i. In the region of large positive gate voltage at zero bias voltage, ej and 
e| are all empty, which means that the two levels are degenerate. Therefore 




Fig. 28 (Color) Nested stability diagram of a parallel DSJ with parameters ei jCT = 
-1.8 eV, £2, CT = -0.3 eV, Ui = 3.6 eV, U 2 = 0.6 eV t = 0.001 eV, 7* = 7^ = 0.04 eV, 
-yl = 7^ = 0.05 eV, V g2 = V gl /2 = V g /2 and Vr = -V h = V blas /2. See discussion in the 
text. 
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adding a bias voltage, first, there will be two channels (ej and ej_) opened 
for current (thick lines). After then, one level e CT (spin-up or spin-down) is 
occupied, while the other obtains a shift for Coulomb interaction: e s — ► e s +U . 
Therefore, when the bias voltage is further increased to make the Fermi- 
window boundary meeting level e s + U , only one channel is opened for the 
current, which results in the weak lines in Fig. 1281 which is the characteristic 
of CB. The similar case appears in the region of large negative gate voltages. 

Finally, we here introduced a powerful Ansatz for the lesser Green func- 
tion, which is consistent with both the Dyson-equation approach and the 
equation-of-motion approach. By using this Ansatz together with the stan- 
dard equation-of-motion technique for the retarded and advanced Green func- 
tions, we obtained the NEGF for both the single and the double site junctions 
in the Coulomb blockade regime at finite voltages and calculated the trans- 
port observables. The method can be applied to describe self-consistently 
transport through single molecules with strong Coulomb interaction and ar- 
bitrary coupling to the leads. 

To test our method, we here analyzed the CB stability diagrams for a SSJ 
and a DSJ. Our results are all consistent with the results of experiments and 
the master-equation approach. We showed, that the improved lesser Green 
function gives better results for weak molcculc-to-contact couplings, where a 
comparison with the master equation approach is possible. 

For the serial configuration of a DSJ, such as a donor/acceptor rectifier, 
the I-V curves maintain a diode-like behavior, as it can be already inferred 
by coherent transport calculations. [265J Besides, wc predict that as a result 
of charging effects, one conductance peak will be split into three peaks and 
one dip, and one current maximum into two. For a DSJ parallel configuration, 
due to different charging energies on the two dot sites, the stability diagrams 
show peculiar nesting characteristics. 
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Jf.,2 Nonequilibrium vibrons 

Though the electron-vibron model described in the Section II has a long 
history, the many questions it implies are not answered up to now. While 
the isolated electron-vibron model can be solved exactly by the so-called po- 
laron or Lang-Firsov transformation [93 [96l [97] , the coupling to the leads 
produces a true many-body problem. The inelastic resonant tunneling of 
single electrons through the localized state coupled to phonons was first con- 
sidered in Refs. [98[ [93 1 1001 HOlj . There, the exact solution in the single- 
particle approximation was derived, ignoring completely the Fermi sea in 
the leads. At strong electron-vibron couplings and weak couplings to the 
leads, satellites of the main resonant peak are formed in the spectral func- 
tion (Fig. rTTj) . The number of the relevant side-bands is determined by the 
well known Huang- Rhys factor [292] g = (X/luq) 2 . The question which re- 
mains is whether these side-bands can be observed in the differential conduc- 
tance, when the coupling to all electrons in the leads should be taken into 
account simultaneously. New theoretical treatments were presented recently 

in Refs. pm ma em qds ens una n® nm ma nsa inn ma eh es 

[841 [T241 [2751 [2761 [T64l [2091 [T25l [2791 [2781 [2771 IT331 [TT2] . 

In parallel, the theory of inelastic scanning tunneling spectroscopy was 
developed [Mil [1621 0331 [IIll [mi E3 USB]- For a recent review of the 
electron-vibron problem and its relation to charge transport at the molecular 
scale see Ref. |164) . Note the related problem of quantum shuttle [1431 11451 
Hg5tH37| . 

Many interesting results by the investigation of quantum transport in the 
strong electron-vibron coupling limit has been achieved with the help of the 
master equation approach [104[ 11061 1107|. 1108] I109j . This method, however, is 
valid only in the limit of very weak molccule-to-lead coupling and neglects all 
spectral effects, which are the most important at finite coupling to the leads. 

4.2.1 Nonequilibrium Dyson-Keldysh method 

(i) The model electron-vibron Hamiltonian 

We use the minimal transport model described in the previous sections. For 
convenience, we present the Hamiltonian here once more. The full Hamilto- 
nian is the sum of the molecular Hamiltonian Hm, the Hamiltonians of the 
leads -Hfl(L), the tunneling Hamiltonian Ht describing the molccule-to-lead 
coupling, the vibron Hamiltonian Hy including electron-vibron interaction 
and coupling of vibrations to the environment (describing dissipation of vi- 
brons) 

H = H M + Hv + H L + H R + H T . (405) 
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A molecule is described by a set of localized states |a) with energies e a 
and inter-orbital overlap integrals t a p by the following model Hamiltonian: 

H ( m ] - E ( £ « + e ^(*)) <k d ° + E toP&dp- (406) 

Vibrations and the electron-vibron coupling are described by the Hamil- 
tonian [ml HM HH [T24] {h = l) 

#V = E "A** + E E A «/?K + 4)4 V (407) 

9 a/3 9 

Here vibrations are considered as localized phonons and q is an index la- 
beling them, not the wave-vector. The first term describes free vibrons with 
the energy tu q . The second term represents the electron-vibron interaction. 
We include both diagonal coupling, which describes a change of the electro- 
static energy with the distance between atoms, and the off-diagonal coupling, 
which describes the dependence of the matrix elements t a p over the distance 
between atoms. 

The Hamiltonians of the right (R) and left (L) leads read 

h 1=L (r) = E( e ^ + e W) c Lr c itoJ ( 408 ) 

(fii(t) are the electrical potentials of the leads. Finally, the tunneling Hamil- 
tonian 

Ht= E E (Wi^+Lc) (409) 

i—L,R ka,a 

describes the hopping between the leads and the molecule. A direct hopping 
between two leads is neglected. 

(ii) Kcldysh-Dyson equations and self-energies 

We use the nonequilibrium Green function (NGF) method, as introduced in 
Section III. The current in the left (i = L) or right (i = R) contact to the 
molecule is described by the expression 

J ^ = l/ ^Tr{r i (e- e ^)(G<(e)+/°(e-e^)[G«(e)-G A (6)])} 1 

. (410) 

where f®(e) is the equilibrium Fermi distribution function with chemical po- 
tential /ii, and the level- width function is 

r i= z,(R)(e) = r ial3 (e) = 2nJ2 v ika,i3V* klja S(e - e ika ). 

fc<7 
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The lesser (retarded, advanced) Green function matrix of a nonequilibrium 
molecule G <( - R - A ' ) = G^ R ' A ^ can be found from the Dyson-Keldysh equations 
in the integral form 

G*(e) = G*(e) + G*(e)£«(e)G«(e), (411) 

G<(e)=G«(e)£<(e)G^( e ), (412) 

or from the corresponding equations in the differential form (see e.g. Refs. (1231 
1124) and references therein). 
Here 

S*'< = £f + s£ <(T) + S *.«v) (413) 

is the total self-energy of the molecule composed of the tunneling (coupling 
to the left and right leads) self-energies 

- ^ (T) = E {y^^<V jk ^ , (414) 



J j=LM 

ha 



R,<(V) 



and the vibronic self-energy £ fl '<( y ) = JJ a 

For the retarded tunneling self-energy one obtains 

Sf (T) (e) = Ai(e - e^) - ^(e - e^), (415) 

where Aj is the real part of the self-energy, which usually can be included in 
the single-particle Hamiltonian H]$ , and Ti describes level broadening due 
to coupling to the leads. For the corresponding lesser function one finds 

£< (T) (e) = il\(e - e^)/?(c - e Vi ). (416) 

In the standard self-consistent Born approximation, using the Keldysh 
technique, one obtains for the vibronic self-energies [1061 11211 11201 11221 11171 

m eii mil 



i k - 

quj 



E*W(e) = iW^(M*Gf_ u M«ZP 

q J 71 

+M"Gf_ u M^« - 2£« =0 M«Tr [G<M?] ) , (417) 
£<W(e)=i£ / ^M*G<_JA*D 



)< 

'qui l 



(418) 



where G K = 2G<+G K -G A is the Keldysh Green function, and M« = M^. 

If vibrons are noninteracting, in equilibrium, and non-dissipative, then the 
vibronic Green functions write: 
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D R (q,u) = 1 —— 1 (419) 

UJ — U!„ + l0 + UJ + UJ„ + l0 + 



D< (q, u)=- 2iri [(/° K) + l)6(u> + uj q ) + f B {w q )5{w - w,)] , (420) 
where the equilibrium Bose distribution function is 

/£M = 7 L 7 - (421) 

exp ) — 1 

In the Migdal model the retarded vibron function is calculated from the 
Dyson-Keldysh equation 

D R M = - 2 2 ^ v (422) 

-<J* - 2uj q n H (q,uj) 

where II(q,uj) is the polarization operator (boson self-energy). The equation 
for the lesser function (quantum kinetic equation in the integral form) is 

{nil - n*)D< - (d r - Dfjn< = o, (423) 

this equation in the stationary case considered here is algebraic in the fre- 
quency domain. 

The polarization operator is the sum of two parts, environmental and 
electronic: II*>< = <J <(cnv) + 77£ <(ol) . 

The environmental equilibrium part of the polarization operator can be 
approximated by the simple expressions 

n R ^\q ) u) = - l - 7q sign(u J ), (424) 

n<™ v \q,tu) = - l79 /£MsignM, (425) 

where j g is the vibronic dissipation rate, and f B (u>) is the equilibrium Bose- 
Einstein distribution function. 

The electronic contribution to the polarization operator within the SCBA 

is 

n R ^(q, u) = -i j |^Tr (M 9 G<M"G^ + M q G R M q Gf_ u ) , (426) 

77< cl > (q, u) = -if ^Tr (M«G e <M«G £ >_ w ) . (427) 

We obtained the full set of equations, which can be used for numerical 
calculations. 
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4.2.2 Single-level model: spectroscopy of vibrons 

The isolated single-level electron-vibron model is described by the Hamilto- 
nian 

H M +v = (eo + e(p )d)d + w a!a + A (a f + a) d) d, (428) 

where the first and the second terms describe the free electron state and the 
free vibron, and the third term is electron-vibron minimal coupling interac- 
tion. 

The electrical potential of the molecule tpo plays an important role in 
transport at finite voltages. It describes the shift of the molecular level by 
the bias voltage, which is divided between the left lead (tip), the right lead 
(substrate), and the molecule as ipo = ipn + r)((fL — ^Pr) .293]. We assume 
the simplest linear dependence of the molecular potential (j] = const), but 
its nonlinear dependence |294j can be easily included in our model. 

Here we assume, that the vibrons are in equilibrium and are not excited 
by the current, so that the self-consistent Born approximation is a good 
starting point. The vibron Green function are assumed to be equilibrium 
with the broadening defined by the external thermal bath, see for details 

Refs. pniiiTnimrii^iT24] . 

For the single-level model all equations are significantly simplified. Com- 
bining Jl and Jr the expression for the current can be written for energy 
independent (wide-band limit) as 

J = Ij^Jjr J deA(e) [f(e - e^ L ) f(e - e<p R )] . (429) 

It looks as simple as the Landauer-Biittiker formula, but it is not trivial, 
because the spectral density A(e) = — 2lmG R (e) now depends on the distri- 
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Fig. 29 (Color online) Spectral function at different electron-vibron couplings: X/luq = 0.4 
(black), X/u>o = 1.2 (blue, dashed), and A/ojo = 2 (red); at to/^O = 5; -^iA^O = ^r/^O 
J. I In the insert the spectral function at X/u>o = 1.2 is shown at finite voltage, when the 
level is partially filled. Energies are in units of Huiq. 
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Voltage 

Fig. 30 Differential conductance of a symmetric junction (rj = 0.5, fjj = 7~£,) at different 
molecule-to-lead coupling, from r^/uiQ = 0.1 (lower curve) to r^/uQ = 10 (upper curve), 
\/u)(i = 1, eo/ujQ = 2. Voltage is in the units of hujo/e. 



bution function of the electrons in the fluctuating molecule and hence the 
applied voltage, ifL = —tpn, = V/2 [123] . Indeed, G R (e) can be found from 

= : rW /uvr A/9 ' ( 43 °) 

e-eo- 2J H ^ v '{e) + i(F L + r R )/2 

where E R( - V \e) is a functional of the electron distribution function inside 
a molecule. Actually, the lesser function G < (e) is used in the quantum ki- 
netic formalism as a distribution function. In the single-level case the usual 
distribution function can be introduced through the relation 

G<(e)=iA(e)f(e). (431) 

Note the essential difference between symmetric (II = r R ) and asymmet- 
ric junctions. It is clear from the nonintcracting solution of the transport 
problem. Neglecting for a moment the vibron self-energies, we obtain the 
nonintcracting distribution function 

= r LfU e ~ ei PL) + r R f R (e - e<p R ) ^ 32 ^ 

-Tl + r R 

For strongly asymmetric junctions (e.g. Ti, <C r R ) the distribution function 
remains close to the equilibrium function in the right lead f R (e — eip R ), thus 
essentially simplifying the solution. While for symmetric junctions the dis- 
tribution function has the double-step form and is very different from the 
equilibrium one. 

A typical example of the spectral function at zero voltage is shown in 
Fig. [29] At finite voltage it should be calculated self-consistently. In the insert 




Voltage 

Fig. 31 Differential conductance of an asymmetric junction [ij = 0, 7jj = 20/Y) at 
different molecule-to-lead coupling, from Fr/loo = 0.2 (lower curve) to -Tfl/oJo = 4 (upper 
curve), A/ojq = 2, eoA^O = 5- The voltage is in the units of hu>o/e 



the spectral function of the symmetric junction at finite voltage is shown, it 
is changed essentially because the distribution function is changed. 

Let us discuss a general picture of the vibronic transport in symmetric 
and asymmetric singlc-molccule junctions, provided in experiments with the 
molecular bridges and STM-to-molcculc junctions, respectively. The differen- 
tial conductance, calculated at different molccule-to-lead coupling, is shown 
in Fig. [30] (symmetric) and Fig. [3T] (asymmetric). At weak coupling, the vi- 
bronic side-band peaks are observed, reproducing the corresponding peaks in 
the spectral function. At strong couplings the broadening of the electronic 
state hides the side-bands, and new features become visible. In the symmet- 
ric junction, a suppression of the conductance at V ~ ±fkjQ takes place as 
a result of inelastic scattering of the coherently transformed from the left 
lead to the right lead electrons. In the asymmetric junction (Fig. I3T|) . the 
usual IETS increasing of the conductance is observed at a negative voltage 
V ~ — hujOi this feature is weak and can be observed only in the incoherent 
tail of the resonant conductance. We conclude, that the vibronic contribu- 
tion to the conductance can be distinguished clearly in both coherent and 
tunneling limits. 

Now let us discuss the particular situation of STS experiments [32] [33] [35] 
I36j . Here we concentrate mainly on the dependence on the tip-to-molcculc dis- 
tance [33J . When the tip (left lead in our notations) is far from the molecule, 
the junction is strongly asymmetric: Jj, <C Tr and r\ — > 0, and the con- 
ductance is similar to that shown in Fig.[3T] When the tip is close to the 
molecule, the junction is approximately symmetric: 1^ w _Tr and r\ 0.5, 
and the conductance curve is of the type shown in Fig.[30j We calculated the 
transformation of the conductance from the asymmetric to symmetric case 
(Fig. [32]) . It is one new feature appeared in asymmetric case due to the fact 
that we started from a finite parameter r\ = 0.2 (in the Fig.[3T]?7 = 0), namely 
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Fig. 32 (Color) Differential conductance at different molcculc-to-STM coupling (sec the 
text), from asymmetric junction with r^/uto = 0.025, -T^/i^o = 0.5 and r\ = 0.2 (lower 
curve, blue thick line) to symmetric junction with J^/ujo = Fr/wo = 0.5 and r\ = 0.5 
(upper curve, red thick line), X/loq = 1> eo/wo = 2. Voltage is in the units of hu)o/e 



a single peak at negative voltages, which is shifted to smaller voltage in the 
symmetric junction. The form and behavior of this peak is in agreement with 
experimental results [55] , 

In conclusion, at weak molecule-to-lead (tip, substrate) coupling the usual 
vibronic side-band peaks in the differential conductance are observed; at 
stronger coupling to the leads (broadening) these peaks are transformed into 
step-like features. A vibronic-induced decreasing of the conductance with 
voltage is observed in high-conductance junctions. The usual IETS feature 
(increasing of the conductance) can be observed only in the case of low off- 
resonant conductance. By changing independently the bias voltage and the 
tip position, it is possible to determine the energy of molecular orbitals and 
the spectrum of molecular vibrations. In the multi-level systems with strong 
electron-electron interaction further effects, such as Coulomb blockade and 
Kondo effect, could dominate over the physics which we address here; these 
effects have to be included in a subsequent step. 

4.2.3 Multi-level model: nonequilibrium vibrons 

Basically there are two main nonequilibrium effects: the electronic spectrum 
modification and excitation of vibrons (quantum vibrations). In the weak 
electron-vibron coupling case the spectrum modification is usually small 
(which is dependent, however, on the vibron dissipation rate, temperature, 
etc.) and the main possible nonequilibrium effect is the excitation of vibrons 
at finite voltages. We have developed an analytical theory for this case |124j . 
This theory is based on the self-consistent Born approximation (SCBA), 
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which allows to take easily into account and calculate nonequilibrium dis- 
tribution functions of electrons and vibrons. 

If the mechanical degrees of freedom are coupled strongly to the envi- 
ronment (dissipative vibron), then the dissipation of molecular vibrations is 
determined by the environment. However, if the coupling of vibrations to the 
leads is weak, we should consider the case when the vibrations arc excited by 
the current flowing through a molecule, and the dissipation of vibrations is 
also determined essentially by the coupling to the electrons. Here , we show 
that the effects of vibron emission and vibronic instability are important 
especially in the case of electron- vibron resonance. 

We simplify the equations and obtain some analytical results in the vi- 
bronic quasiparticle approximation, which assumes weak electron-vibron cou- 
pling limit and weak external dissipation of vibrons: 

7g * = lq - 2lmn R (Lu q ) < uj q . (433) 

So that the spectral function of vibrons can be approximated by the Dirac <5, 
and the lesser function reads 

D< (q, u) = -2ni [(N q + l)6(u> + w,) + N q S(uj - u> q )] , (434) 

where N q is (nonequilibrium) number of vibrations in the g-th mode. So, 
in this approximation the spectrum modification of vibrons is not taken 
into account, but the possible excitation of vibrations is described by the 
nonequilibrium N q . The dissipation of vibrons is neglected in the spectral 
function, but is taken into account later in the kinetic equation for N q . 
A similar approach to the single-level problem was considered recently in 
[TT31 fTT4l Ell niSl HQ61 Ell EH]. The more general case with broadened 
equilibrium vibron spectral function seems to be not very interesting, be- 
cause in this case vibrons are not excited. Nevertheless, in the numerical 
calculation it can be easy taken into consideration. 

From the general quantum kinetic equation for vibrons, we obtain in this 
limit 

This expression describes the number of vibrons N q in a nonequilibrium 
state, N® = f%(w q ) is the equilibrium number of vibrons. In the linear ap- 
proximation the polarization operator is independent of N q and —2lmII R (u> q ) 
describes additional dissipation. Note that in equilibrium N q = 2V? because 
lmll < (u! q ) = 2lmII R (uj q )f t g(uj q ). See also detailed discussion of vibron emis- 
sion and absorption rates in Rcfs. |113|4114| I115|, [TT6] . 

For weak electron-vibron coupling the number of vibrons is close to equilib- 
rium and is changed because of vibron emission by nonequilibrium electrons, 
N q is roughly proportional to the number of such electrons, and the distri- 
bution function of nonequilibrium electrons is not change essentially by the 
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interaction with vibrons (perturbation theory can be used). The situation 
changes, however, if noncquilibrium dissipation — 2lmII R (ui q ) is negative. In 
this case the number of vibrons can be essentially larger than in the equi- 
librium case (vibronic instability), and the change of electron distribution 
function should be taken into account self-consistcntly. 
In the stationary state the nonlinear dissipation rate 



is positive, but the nonequilibrium contribution to dissipation — 2lmII R (u} q ) 
remains negative. 

Additionally to the vibronic quasiparticle approximation, the electronic 
quasiparticle approximation can be used when the coupling to the leads is 
weak. In this case the lesser function can be parameterized through the num- 
ber of electrons F n in the eigenstates of the noninteracting molecular Hamil- 
tonian Hw 



we introduce the unitary matrix S, which transfer the Hamiltonian H = 
Hjrf a p hito the diagonal form H = S _1 HS, so that the spectral function of 
this diagonal Hamiltonian is 



7g * = lq -2lmn R (w q ) 



(436) 




(437) 



A Sri (e) 



27T(5(e - es)5g v , 



(438) 



where es are the eigenenergics. 
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Fig. 33 (Color) Vibronic emission in the symmetric multilevel model: voltage-current 
curve, differential conductance, and the number of excited vibrons in the off-resonant 
(triangles) and resonant (crosses) cases (details see in the text). 
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Note that in the calculation of the self-energies and polarization operators 
we can not use ^-approximation for the spectral function (this is too rough 
and results in the absence of interaction out of the exact electron-vibron 
resonance). So that in the calculation we use actually ()437|) with broad- 
ened equilibrium spectral function. This approximation can be systematically 
improved by including noncquilibrium corrections to the spectral function, 
which are important near the resonance. It is important to comment that 
for stronger electron-vibron coupling vibronic side-bands are observed in the 
spectral function and voltage-current curves at energies e,5 ± nu) q , we do not 
consider these effects in the rest of our paper and concentrate on resonance 
effects. 

After correspondingly calculations we obtain finally 



7^°-£^M" 9 )W-i) 



(439) 



where coefficients K v g are determined by the spectral function and electron- 
vibron coupling in the diagonal representation 



de - 



2vr 



M^Assie-LO^M^ie), 



(440) 



F„ 







Q q F S N q +tts q m- 




^7777+^7777 + Ylqr) 






-N q ) 



(441) 



(442) 



here T irjn and are the level width matrix in the diagonal representation 
and Fermi function at energy e r; — eipt. 

These kinetic equations are similar to the usual golden rule equations, but 
are more general. 

Now let us consider several examples of vibron emission and vibronic in- 
stability. 

First we consider the most simple case, when the instability is not possible 
and only vibron emission takes place. This corresponds to a negative imagi- 
nary part of the electronic polarization operator: lraJI R ' (cjq) < 0. From the 
Eq. (|440p one can see that for any two levels with the energies e n > the 
coefficient k v s is larger than Ks-q, because the spectral function ^^(e) has 
a maximum at e = e$. The contribution of K ri s(uo q ){F rj — F$) is negative if 
F n < F$. This takes place in equilibrium, and in noncquilibrium for trans- 
port through symmetric molecules, when higher energy levels are populated 
after lower levels. The example of such a system is shown in Fig.[33J Here 
we consider a simple three-level system (ei = 1, ?2 = 2, 63 = 3) coupled 
symmetrically to the leads = = 0.01). The current-voltage curve 

is the same with and without vibrations in the case of symmetrical coupling 
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Voltage 



Fig. 34 (Color) Vibronic instability in an asymmetric multilevel model: voltage-current 
curve, differential conductance, and the number of excited vibrons (crosses). Dashed line 
show the voltage-current curve without vibrons (details see in the text). 




- 1 *r .° 1 

Voltage 

Fig. 35 (Color) Floating level resonance: voltage-current curve and the number of excited 
vibrons (crosses). Dashed line show the voltage-current curve without vibrons (details see 
in the text). 



to the leads and in the weak electron-vibron coupling limit (if we neglect 
change of the spectral function). The figure shows how vibrons are excited, 
the number of vibrons Ny in the mode with frequency ujq is presented in two 
cases. In the off-resonant case (green triangles) Ny is very small comparing 
with the resonant case (ujq = £2 — ?i, red crosses, the vertical scale is changed 
for the off-resonant points). In fact, if the number of vibrons is very large, the 
spectral function and voltage-current curve are changed. We shall consider 
this in a separate publication. 

Now let us consider the situation when the imaginary part of the electronic 
polarization operator can be positive: lmII R (uj q ) > 0. Above we considered 
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the normal case when the population of higher energy levels is smaller than 
lower levels. The opposite case F 2 > F\ is known as inversion in laser physics. 
Such a state is unstable if the total dissipation 7* (|436j) is negative, which is 
possible only in the nonstationary case. As a result of the instability, a large 
number of vibrons is excited, and in the stationary state 7* is positive. This 
effect can be observed for transport through asymmetric molecules, when 
higher energy levels are populated before lower levels. The example of a such 
system is shown in Fig. [34] It is the same three-level system as before, but 
the hrst and second levels are coupled not symmetrically to the leads (Ili = 
0.001, r R1 = 0.1, r L2 = O.l, r R2 = O.OOl). The vibron couple resonantly 
these levels (uj q = £2 — ei). The result is qualitatively different from the 
symmetrical case. The voltage-current curve is now asymmetric, a large step 
corresponds to the resonant level with inverted population. 

Note the importance of the off-diagonal electron-vibron coupling for the 
resonant effects. If the matrix M in the cigen-state representation is diagonal, 
there is no resonant coupling between different electronic states. 

Finally, let us consider the important case, when initially symmetric 
molecule becomes asymmetric when the external voltage is applied. The rea- 
son for such asymmetry is simply that in the external electric field left and 
right atoms feel different electrical potentials and the position of the levels 
e« = Cq '' + ei/3 Q is changed (float) with the external voltage. The example of 
a such system is shown in Fig. 023 Here we consider a two-level system, one 
level is coupled electrostatically to the left lead 1\ cx ip^, the other level to the 
right lead e 2 oc tp R , the tunneling coupling to the leads also is not symmetri- 
cal (r L1 = 0.1, r R1 = 0.001, r L2 = O.OOl, r R2 = O.l). The frequency of the 
vibration, coupling these two states, is ujq = 1. When we sweep the voltage, 
a peak in the voltage-current curve is observed when the energy difference 
1\ — e 2 tx eV is going through the resonance £\ — e 2 loq. 
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4-3 Coupling to a vibrational contiuum: dissipation 
and renormalization 

4.3.1 The model Hamiltonian 

In the previous section we have dealt with a simple, but nevertheless phys- 
ically rich, model describing the interaction of an electronic level with some 
specific vibrational mode confined to the quantum dot. We have seen how 
to apply in this case the Keldysh non-equilibrium techniques described in 
Section III within the self-consistent Born and Migdal approximations. The 
latter are however appropriate for the weak coupling limit to the vibrational 
degrees of freedom. In the opposite case of strong coupling, different tech- 
niques must be applied. For equilibrium problems, unitary transformations 
combined with variational approaches can be used, in non-equilibrium only 
recently some attempts were made to deal with the problem. |139j 

In this section we will consider the case of a multi-level electronic system 
in interaction with a bosonic bath (2881 1289] . We will use unitary transforma- 
tion techniques to deal with the problem, but will only focus on the low-bias 
transport, so that strong non-equilibrium effects can be disregarded. Our in- 
terest is to explore how the qualitative low-energy properties of the electronic 
system are modified by the interaction with the bosonic bath. We will see that 
the existence of a continuum of vibrational excitations (up to some cut-off fre- 
quency) dramatically changes the analytic properties of the electronic Green 
function and may lead in some limiting cases to a qualitative modification 
of the low-energy electronic spectrum. As a result, the I-V characteristics at 
low bias may display "metallic" behavior (finite current) even if the isolated 
electronic system does exhibit a band gap. The model to be discussed below 
has been motivated by the very exciting electrical transport measurements 
on short poly(dG)-poly(dC) DNA molecular wires carried out at the group 
of N. Tao some time ago [60]. Peculiar in these experiments was the large 
measured currents -up to 150 nA at 0.8 V- at low voltages, which stood in 
strong contrast to the usually accepted view that DNA should behave as an 
insulator at low applied bias. Further, a power-law length scaling of the linear 
conductance with increasing wire length was demonstrated, indicating that 
long-range charge transport was possible. Since the experiments were carried 
out in an aqueous solution, the possibility of a solvent-induced modification 
of the low-energy transport properties of the wire lied at hand, although 
additional factors like internal vibrations could also play a role. 

The proposed model is based on an earlier work [286j and assumes, within 
a minimal tight-binding picture, that the DNA electronic states can be quali- 
tatively classified into extended (conducting) and localized (non-conducting) 
states. The former may correspond e.g. to the 7r-orbital stack of the base 
pairs, the latter to energetically deeper lying (w.r.t. the frontier orbitals) 
base-pair states or sugar-phosphate backbone states. A further assumption is 
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Fig. 36 (Color) Schematic drawing of a DNA molecular wire in contact with a dissipative 
environment. The central chain (extended states) with TV sites is connected to semiinfinitc 
left (L) and right (R) electronic reservoirs. The bath only interacts with the side chain 
sites (c), which we call for simplicity backbone sites, but which collectively stay for non- 
conducting, localized electronic states. The Hamiltonian associated with this model is given 
by Eqs. 14431 1. 14441 . and 14451 1 in the main text. 

that any modification of the conducting states through the environment only 
takes place through a coupling to the non-conducting set. The tight-binding 
electronic Hamiltonian for N sites can then be written as (see also Fig. 131)1) : 



Hereby Tic and Tib are the Hamiltonians of the extended and localized 
states (called in what follows "backbone" states for simplicity), respectively, 
and Hc-b is the coupling between them. t\\ and t± are hopping integrals 
along the central chain (extended states) and between the localized states 
and the central chain, respectively. If not stated otherwise, the on-site ener- 
gies will be later set equal to zero to simplify the calculations. Notice that 
this model displays a gap in the electronic spectrum roughly proportional to 
the transversal coupling t±. This can be easely seen by looking at the limit 
N —> oo which leads to a periodic system. In this case, the Hamiltonian can 
be analytically diagonalized and two energy dispersion curves are obtained, 



which are given by E±(k) = t|| cos(fc) ± Jt\ +tnCOS 2 (fc). The direct gap 



between the two bands is simply 5 = 2y£^ "^fl" Since this model further 
shows electron-hole symmetry, two electronic manifolds (bands in the limit 
of N — > oo ) containing N states each, are symmetrically situated around the 
Fermi level, which is taken as the zero of energy. 

The gap is obviously temperature independent and furthermore it is ex- 
pected that transport at energies E < S will be strongly suppressed due to 




(443) 
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the absence of electronic states to suport charge propagation. As a result, 
the linear conductance should display a strong exponential dependence as a 
function of the chain length N. In view of this behavior, an immediate issue 
that arises is how stable this electronic structure, i.e. two electronic manifolds 
separated by a gap, is against the coupling to an environment. This is an issue 
which reaches farther than the problem of charge transport in DNA wires, 
since it addresses the interaction of an open quantum mechanical system with 
a countable number of electronic energy levels to a continuum of states ( "uni- 
verse" ) . A generic example of such a situation is the measurement process in 
quantum mechanics. It is well-known that the interaction with complex envi- 
ronments is a source of dissipation and decoherence in quantum mechanical 
systems. [71] Concerning more specifically the case of DNA (and proteins), 
there is broad experimental evidence that the molecule dynamics follows the 
solvent dynamics over a broad temperature range. Especially, conformational 
changes, low-energy vibrational excitations and the corresponding tempera- 
ture dependences turn out to be very sensitive to the solvents dynamics. |295j 
We will thus consider the vibrational degrees of freedom of counterions and 
hydration shells of the solvent as a dynamical bath able to break the elec- 
tronic phase coherence and additionally to act as a dissipative environment. 
We do not consider specific features of the environment but represent it in a 
generic way by a bosonic bath of M harmonic oscillators. Then, the previous 
Hamiltonian can be extended to: 

H w = n e i + J2^cBiB a + J2^a4cj(B a + Bl) = n ei + n B +n c . B , (444) 

where 7Yb and 7Y c _b are the phonon bath Hamiltonian and the (localized) 
state-bath interaction, respectively. B a is a bath phonon operator and \ a 
denotes the electron-phonon coupling. Note that we assume a local coupling 
of the bath modes to the electronic density at the side chain. Later on, the 
thermodynamic limit (M — ► oo) in the bath degrees of freedom will be carried 
out and the corresponding bath spectral density introduced, so that at this 
stage we do not need to further specify the set of bath frequencies Q a and 
coupling constants A Q . Obviously, the bath can be assumed to be in thermal 
equilibrium and be described by a canonical partition function. 

To complete the formulation of the model, we have to include the inter- 
action of the system with electronic reservoirs in order to describe charge 
transport along the same lines as before. We assume, as usual, a tunnel-type 
Hamiltonian with the form: 

H=H W + e ^ d la d ^ + ( VkA ^ka fe l + H - C 

keL,R,er keL,cr 

+ ]T (V k , N 4 a b N + H.c.) =H W + H L/R + Wl-c + Hr-c- (445) 

kefl.cr 
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The Hamiltonian of Eq. (|445p is the starting point of our investigation. 
For a weak charge-bath coupling, a perturbative approach similar to the 
second order Born approximation, as described in the previous section can be 
applied. We expect, however, qualitative new effects rather in the opposite 
limit of strong coupling to the bath. To deal with this problem, a unitary 
transformation, the Lang-Firsov (LF) transformation, can be performed on 
the Hamiltonian of Eq. (|445j) . which allows to eliminate the linear charge- 
vibron interaction 7i c -B- In the limiting case of an isolated system with a 
single electron (or hole) this transformation becomes exact and allows for a 
full decoupling of electronic and vibronic propagators, see e.g. Ref. (37]. In 
the present case, this transformation is not exact and further approximations 
have to be introduced in order to make the problem tractable. 

The generator of the LF transformation is given by 

S = Y J ^ a /n a )c]c 1 {B a -Bi) 

and = —S. In the transformed Hamiltonian H = c s Hc~ s the linear 
coupling to the bath is eliminated. One should notice that in H only the 
"backbone" part of the Hamiltonian is modified since the conducting state 
operators bi as well as the lead operators are invariant with respect to 
the above transformation. The new Hamiltonian reads: 

H = H C + Wl/r + H b + H l/r . c + (e-A)J2 4 C J ~ *-L [ b h X + Hc - 



3 



X = cxp 



As a result of the LF we get a shift of the onsite energies (polaron shift 
or reorganization energy in electron transfer theory) and a renormalization 
of both the tunneling and of the transversal coupling Hamiltonian via the 
bosonic operators X. There is also an additional electron-electron interaction 
term which we will not be concerned with in the remaining of this section and 
is thus omitted. Since we are mainly interested in qualitative statements, we 
will assume the wide-band approximation in the coupling to the electrodes 
which is equivalent to substituting the electrode self-energies by a purely 
imaginary constant, i.e. £l,r ~ — i ^l,r- We are thus not interested in specific 
features of the electrode electronic structure. 

To further proceed, let us now introduce two kinds of retarded thermal 
Green functions related to the central chain Gje(t) and to the "backbones" 
Pje{t), respectively (taking h = 1): 
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G je (t,t') = -iG(t-t')( b 3 {t),b\{t') 



(447) 



Pjt(t,t') = -iG(t-t'){ C] (t)X(t),cl(t')X\t') 



where is the Heaviside function. Notice that the P-Grccn function doe not 
have a pure electronic character but also contains the bath operators X. For 
a full out-of-cquilibrium calculation, the full Kcldysh formalism including 
lesser- and greater-GF would also be needed. However, as we will briefly 
show below, the final expression for the electrical current at low applied 
voltages and for small transversal coupling t± will only include the retarded 
propagators. 

We now use the equation of motion technique (EOM) to obtain an expres- 
sion for the GF Gji(t). We first remark that in the time domain two EOM 
can be written, depending on which time argument in the double-time GF 
the time derivative will act. One thus obtains in general: 



id t G(t, t') = ( [b(t),tf(t')] + )6(t f) + (([b(t), H] |&t(f))). 
G(t, 0[-i d t ,\ = ( [b(t), &t(f)] + ) S (t t') {{b{t)\ [&t(f), H ] )). 

The EOM for the GF Gji(t) reads then in the energy space: 

Y, [G^{E)] tn G nj {E) = ^ - t±((ciX\b])) 



(448) 



hi 



>L(R) 



l^k,l(JV)l 



E 

k£L(_R) 



£-e k + iO+ 



-i-fl,R 



In the next step, EOM for the "right" time argument t' of the GF 
Z^{t,t'){{cg{t)X{t)\b](t'))) can be written. This leads to: 

ZL(E) [Go '(E)}^ = -t ± ((ctX\4x*)) = -t ± P ej (E) (449) 

tn 

Inserting Eq. (|449p into Eq. (|448p we arrive at the matrix equation: 

G(E) = G (E) + G Q (E)-£ B (E)G Q (E), 

which can be transformed into a Dyson-like equation when introducing the 
irreducible part E B (£) = ^W(E) + £jf (£)G (£)£g r (£) + ...: 



G(E) = G (E) + Go(E)-££(E)G(E), 



(450) 
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or cquivalcntly: 

G- 1 (E) = G- 1 (E)-tlP(E) (451) 
Go \E) =E1-H C - £l(E) - S R (E). 

Sg r (£) = t\T(E) is the crucial contribution to the GF since it contains 
the influence of the bosonic bath. Note that Sg r (£') includes the transversal 
hopping t± to all orders, the leading one being t\. 

In the next step, an expression for the electrical current flowing through 
the system must be derived. Using the results of Sec. 2, we can directly write 
the following expression: 

I=^f dE Tr(/ L (E) - f R (E)) t(E) 

+ ti^ J dE {Tr[S> P< - £< P>] -(L~R)}. (452) 

The first summand has the same form as Landauer's expression for the cur- 
rent with an effective transmission function t(E) = Tr[G+r R Gr R ]. However, 
the reader should keep in mind that the GFs appearing in this expression do 
contain the full dressing by the bosonic bath and hence, t(E) does not de- 
scribe elastic transport. The remaining terms contain explicitly contributions 
from the bath. It can be shown after some transformations that the leading 
term is proportional to (t^_) 2 so that within a perturbative approach in t± 
and at low bias it can be approximately neglected. We therefore remian with 
the exression I = ^ JdETr(f L (E) - f R {E)) t(E) to obtain the current. 

To remain consistent with this approximation, the bath selfenergy should 
also be treated to order more explicitly: 

P ej (t,t') = ((c e (t)X(t)\c](t')X^ (f))) 

« -i9(t {(«(t)ct(f )) <*(t)*t(f)) + (c]{t')ci{t)) (*t(f )*(f)>} 
« -i6 tj 6(t t>) {( Cj -(t)cJ(f)) <Af(*)Aft(f)> + {c](t') Cj (t)) (x\t')X{t))} 
= -i S tj 9{t - t')e - 1 ^ f { (1 - / c )c + f c e } . (453) 

In the previous expression we have replaced the full averages of the "back- 
bone" operators by their zero order values (free propagators), e - *^ = 
(^X{t)X^ (0)) B is a dynamical bath correlation function to be specified later 
on. The average (-) B is performed over the bath degrees of freedom. f c is the 
Fermi function at the backbone sites. In what follows we consider the case of 
empty sites by setting f c = 0. The Fourier transform Pgj{E) reads then: 
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P ej (E)=-iS ej dte i(E+iO+)t e -i{e-A)t r (1 _ /c)e -*(i) + /ce -*(-t) 



(454) 

In order to get closed expressions for the bath thermal averages it is ap- 
propriate to introduce a bath spectral density [71] defined by : 



(455) 



where uj c is a cut-off frequency related to the bath memory time r c ~ w" 1 . 
It is easy to show that the limit lu c — > oo corresponds to a Markovian bath, 
i.e. J(t) ~ JoS(t). Using this Ansatz, <P(t) can be written as: 



dw 



I — e 



-i 



I — cos wt 



(456) 



Although the integral can be performed analytically [7T], we consider <P(t) in 
some limiting cases where it is easier to work directly with Eq. (|456p . 



4.3.2 Limiting cases 

We use now the results of the foregoing section to discuss the electronic 
transport properties of our model in some limiting cases for which analytic 
expressions can be derived. We will discuss the mean-field approximation 
and the weak-coupling regime in the electron-bath interaction as well as to 
elaborate on the strong-coupling limit. Farther, the cases of ohmic (s = I) 
and superohmic (s = 3) spectral densities are treated. 



(i) Mean-field approximation 



The mean-field approximation is the simplest one and neglects bath fluc- 
tuations contained in P{E). The MFA can be introduced by writing the 
phonon operator X as (X) B + 5X in H c . c in Eq. (f446|) . i.e. = 

— 1± J^j "jCj (X) B + H.c. + 0(5X). As a result a real, static and tempera- 
ture dependent term in Eq. (|45I [) is found: 

G- 1 (E) = G Q 1 (E)-tl 



E-e 



i0+ 



1, 



(457) 



where |(A') B | 2 = e 2k ( t ) and k(T) is given by: 



k{T) 



—~J{u)) coth 



2fc B T 



(458) 
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Fig. 37 (Color) Electronic transmission and corresponding current in the mean-field ap- 
proximation for two different temperatures. Parameters: N = 20, Jo/^c = 0.12,tj_/t| = 
0-5,r L/R /t|| =0.5. 



The effect of the MF term is thus to scale the bare transversal hopping t± 
by the exponential temperature dependent factor e~ K ( T \ 

In the case of an ohmic bath, s = 1, the integrand in k(T) scales as 
l/w p , p = 1, 2 and has thus a logarithmic divergence at the lower integration 
limit. Thus, the MF contribution would vanish. In other words, no gap would 
exist on this approximation level. 

In the superohmic case (s = 3) all integrals are regular. One obtains 
A = J duj uj^ 1 J(oj) = r(s—l)Jo = 2 Jo, with -T(s) being the Gamma function 
and k(T) reads: 



k(T) = 

UJ C 



(459) 



Cn{s, z) = X)^°=o( n + i s the Hurwitz ^-function, a generalization of the 
Riemann ^-function. |296j 

It follows from Eq. (17) that k(T) behaves like a constant for low tem- 
peratures {k^T/uj c < 1), k(T) ~ Jq/uj c , while it scales linear with T in the 
high-tcmpcraturc limit [k^T/uj c > 1), n(T) ~ Jq/u> c (1 + 2ksT / lo c )) . 

For Jo ^ and at zero temperature the hopping integral is roughly reduced 
_i° 

to t±e "a which is similar to the renormalization of the hopping in Holstein's 
polaron model [297], though here it is t± rather than tu the term that is 
rescaled. At high temperatures t± is further reduced (n(T) ~ T) so that 
the gap in the electronic spectrum finally collapses and the system becomes 
"metallic" , see Fig. [37l An appreciable temperature dependence can only be 
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Fig. 38 (Color) Electronic transmission and corresponding current in the weak-coupling 
limit with ohmic dissipation (s = 1) in the bath. Parameters: N = 20, Jo/uj c = 0.2, tx/t\\ = 

o.6,r L/R /t N =0.5 



observed in the limit Jq/uj c < 1; otherwise the gap would collapse already 
at zero temperature due to the exponential dependence on Jo. We further 
remark that the MFA is only valid in this regime (Jo/uj c < 1), since for 
Jq/ujc 3> 1 multiphonon processes in the bath, which are not considered in the 
MFA, become increasingly relevant and thus a neglection of bath fluctuations 
is not possible. 



(ii) Beyond MF: weak-coupling limit 

As a first step beyond the mean-field approach let's first consider the weak- 
coupling limit in P(-E'). For Jq/uj c < 1 and not too high temperatures 
(fcBT/w c < f) the main contribution to the integral in Eq. (|456|) comes from 
long times t 3> ui^ 1 ■ With the change of variables z = ujt, <P(t) can be written 
as: 



(460) 




As far as u c t (3uj c this can be simplified to: 
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/* oo 



/ ; „ Boj c 1 — cos z \ , 
x l- C - 12 + 2 — . 461 

Since in the long-time limit the low-frequency bath modes are giving the 
most important contribution we may expect some qualitative differences in 
the ohmic and superohmic regimes. For s = 1 we obtain <P(t) ~ 7r^^j-^(oj c i) 
which leads to (A(s = 1) = Jo): 

G- 1 ^) = Go- 1 ^) - 1\ 1 1, (462) 

S + Jo + itt^/cbT 

i.e. there is only a pure imaginary contribution from the bath. For the sim- 
ple case of N = 1 (a two-states model) one can easily see that the gap 
approximately scales as \fk~oT; thus it grows with increasing temperature. 
This is shown in Fig. [38l where we also see that the intensity of the trans- 
mission resonances strongly goes down with increasing temperature. The gap 
enhancement is induced by the suppression of the transmission peaks of the 
frontier orbitals, i. e. those closest to the Fermi energy. 

For s = 3 and k^T '/tu c < 1, <P(t) takes a nearly temperature independent 
value proportional to Jo/uj c . As a result the gap is slightly reduced (t± — > 
t±e~ J °/ Mc ) but, because of the weak-coupling condition, the effect is rather 
small. 

From this discussion we can conclude that in the weak-coupling limit ohmic 
dissipation in the bath induces an enhancement of the electronic gap while 
superohmic dissipation does not appreciably affect it. In the high-temperature 
limit k-g,T/uj c > 1 a short-time expansion can be performed which yields 
similar results to those of the strong-coupling limit (see next section) , |298j 
so that we do not need to discuss them here. Note farther that the gap 
obtained in the weak-coupling limit is an "intrinsic" property of the electronic 
system; it is only quantitatively modified by the interaction with the bath 
degrees of freedom. We thus trivially expect a strong exponential dependence 
of t{E = E-p), typical of virtual tunneling through a gap. Indeed, we find 
t(E = E F ) ~ cxp(-/3L) with/3-2-31" 1 . 



(iii) Beyond MF: strong coupling limit (SCL) 

In this section we elaborate on the strong-coupling regime, as defined by the 
condition Jq/lu c > 1. In the SCL the main contribution to the time integral 
in Eq. (|456[) arises from short times. Hence a short-time expansion of <P(t) 
may already give reasonable results and it allows, additionally, to find an 
analytical expression for P(-E'). At t <C w" 1 we find, 
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Fig. 39 (Color) Temperature dependence of the real and imaginary parts of P(E) for 
N = 20, Jq/u>c = 10,ixA|| = 0-4, ^l/r/*| | = 0-5- With increasing temperature the slope 
of the real part near E = decreases and the imaginary part broadens and loses intensity. 
A similar qualitative dependence on Jq was found (not shown). 



<P(t) &iAt + (uj c t) 2 k (T) 



(463) 



P tj (E) = -iS l 



dt e i (B-e+i0+)t e -(w c t) 2 Ko(T) 



1 



exp 
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4w> (T) 
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Before presenting the results for the electronic transmission, it is useful 
to first consider the dependence of the real and imaginary parts of P(E) 
on temperature and on the reduced coupling constant Jq/lu c . Both functions 
are shown in Fig. [3^1 We see that around the Fermi level at E = the 
real part is approximately linear, Re P(E) ~ E while the imaginary part 
shows a Lorentzian-like behavior. The imaginary part loses intensity and 
becomes broadened with increasing temperature or Jo, while the slope in 
the real part decreases when k^T or Jo are increased. If we neglect for the 
moment the imaginary part (the dissipative influence of the bath), we can 
understand the consequences of the real part being nonzero around the Fermi 
energy, i.e. in the gap region of the model of Ref. [286j . The solutions of the 
non-linear equation det|(i? — ij_Re P(E))1 — Hc\ = give the new poles 



Green function techniaucs in the treatment of auantum transDort 



117 




Fig. 40 (Color) Upper panel: t(E) with ImP(B) = 0; the intensity of the resonances on 
the central narrow band is strongly dependent on Jq/w c and kgT (not shown). Temper- 
ature dependence of t(E) with full inclusion of P(E) (middle panel) and corresponding 



current (lower panel) for N 
increases with temperature. 



20, J /u c = 5,t_i_/t|| = 0.5,r L/R /t|| = 0.2. The pseudo-gap 



of the Green function of the system in presence of the phonon bath. For 
comparison, the equation determining the eigcnstates without the bath is 
simply det\(E-t\/E)l-7ic\ = 0. It is just the l/E dependence near E = 
that induces the appearance of two electronic bands of states separated by a 
gap. In our present study, however, KeP(E — > 0) has no singular behavior 
and additional poles of the Green function may be expected to appear in the 
low-energy sector. This is indeed the case, as shown in Fig. [40] (upper panel). 
We find a third band of states around the Fermi energy, which we may call 
a polaronic band because it results from the strong interaction between an 
electron and the bath modes. The intensity of this band as well as its band 
width strongly depend on temperature and on Jo. When feeT (or Jo) become 
large enough, these states spread out and eventually merge with the two other 
side bands. This would result in a transmission spectrum similar of a gapless 
system. 

This picture is nevertheless not complete since the imaginary component 
of P(E) has been neglected. Its inclusion leads to a dramatic modification 
of the spectrum, as shown in Fig. [40] (middle panel). We now only see two 
bands separated by a gap which basically resembles the semiconducting-type 
behavior of the original model. The origin of this gap or rather pseudo-gap 
(see below) is however quite different. It turns out that the imaginary part 
of P{E), being peaked around E = 0, strongly suppresses the transmission 
resonances belonging to the third band. Additionally, the frontier orbitals 
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Fig. 41 (Color) Upper panel: Arrhenius plot for t(Ep). Parameters: N = 20, i| = 
0.6eV, <x/i|| = 0.2, i~L/R,/t|| = 0.3. Middle and lower panels: Length dependence of t(Ep) 
at different temperatures for two different strengths of the electron-bath coupling Jo/co c - 
The electronic coupling parameters are the same as in the upper panel. 



on the side bands, i.e. orbitals closest to the gap region, are also strongly- 
damped, this effect becoming stronger with increasing temperature (ImP(£) 
broadens). This latter effect has some similarities with the previously dis- 
cussed weak-coupling regime. Note, however, that the new electronic mani- 
fold around the Fermi energy does not appear in the weak-coupling regime. 
We further stress that the density of states around the Fermi level is not ex- 
actly zero (hence the term pseudo-gap); the states on the polaronic manifold, 
although strongly damped, contribute nevertheless with a finite temperature 
dependent background to the transmission. As a result, with increasing tem- 
perature, a crossover from "semiconducting" to "metallic" behavior in the 
low-voltage region of the I-V characteristics takes place, see Fig. 00] (lower 
panel). The slope in the I-V plot becomes larger when t± is reduced, since 
the side bands approach each other and the effect of Im P(E) is reinforced. 

In Fig. 22 (top panel) an Arrhenius plot of the transmission at the Fermi 
energy is shown, which suggests that activated transport is the physical mech- 
anism for propagation at low energies. Increasing the coupling to the phonon 
bath makes the suppression of the polaronic band around E = less effective 
(lmP(E ~ 0) decreases) so that the density of states around this energy be- 
comes larger. Hence the absolute value of the transmission will also increase. 
On the other side, increasing t± leads to a reduction of the transmission at 
the Fermi level, since the energetic separation of the side bands increases 
with t±. 
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A controversial issue in transport through molecular wires is the actual 
length dependence of the electron transfer rates or correspondingly, of the 
linear conductance. This is specially critical in the case of DNA nanowires 
[59l [56l E2] . Different functional dependences have been found in charge trans- 
fer experiments ranging from strong exponential behavior related to superex- 
change mediated electron transfer [56] to algebraic dependences typical of 
thermal activated hopping |59( I57| . As far as transport experiments are con- 
cerned, the previously mentioned experiments at the group of N. Tao [60j 
reported an algebraic length dependence of the conductance for poly(GC) 
oligomers in solution. We have investigated the length dependence of t(Ep) 
and found for the strong dissipative regime Jq/lu c > 1, an exponential law for 
energies close to Ep, t(E-p) ~ exp(— 7-ZV). At the first sight, this might be not 
surprising since a gap in the spectrum does exist. Indeed, in the absence of 
the bath, i.e. with an intrinsic gap, we get decay lengths 7 co h of the order of 
2 A . However, as soon as the interaction with the bath is included, we find 
values of 7 much smaller than expected for virtual tunneling, ranging from 
0.15 A to 0.4 A . Additionally, 7 is strongly dependent on the strength 
of the electron-bath coupling Jq/uj c as well as on temperature; 7 is reduced 
when J$/lu c or k^T increases, since in both cases the density of states within 
the pseudo-gap increases. Remarkably, a further increase of the electron-bath 
coupling eventually leads to an algebraic length dependence, see lower panel 
of Fig. HU 

The studies presented in this section indicate that the presence of a com- 
plex environment, which induces decoherence and dissipation, can dramat- 
ically modify the electronic response of a nanowire coupled to electrodes. 
Electron transport on the low-energy sector of the transmission spectrum is 
supported by the formation of (virtual) polaronic states. Though strongly 
damped, these states manifest nonetheless with a finite density of states in- 
side the bandgap and mediate thermally activated transport. 
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5 Conclusions and Perspectives 

In this chapter we have reviewed the method of noncquilibrium Green func- 
tions and few selected applications to problems related with charge transport 
at the molecular scale. Hereby we have only focused on minimal model Hamil- 
tonian formulations which build a very appropriate starting point to illustrate 
the power and range of validity of such techniques. We have showed how this 
approach can be used to deal with a variety of physical systems, covering 
both nonintcracting and interacting cases. Thus, so different issues as co- 
herent transport, Coulomb blockade phenomena, charge- vibron interaction, 
coupling to dissipative environments, and the Kondo effect (not addressed 
in this review) can be in principle treated on the same footing. Specially, 
the existence of well-developed diagrammatic techniques allows for a system- 
atic treatment of interactions in nanoscale quantum systems. For the sake of 
space, we did not deal with applications of NGF techniques to spin-dependent 
transport, a field that has been increasingly attracting the attention of the 
physical community in the past years due to its potential applications in 
quantum information theory and quantum computation 299, 300 . For the 
same reason, the implementation of NGF into first-principle based approaches 
was not discussed neither. This is nevertheless a crucial methodological issue, 
since system-specific and realistic information about molecule-metal contact 
details, charge transfer effects, modifications of the molecular electronic struc- 
ture and configuration upon contacting, the electrostatic potential distribu- 
tion in a device, etc can only be obtained via a full ab initio description of 
transport. For charge transport through noninteracting systems this has been 
accomplished some years ago by combining NGF with DFT methods. The in- 
clusion of interactions, however, represents a much stronger challenge and has 
been mainly carried out, within the self-consistent Born-approximation, for 
the case of tunneling charges coupling to vibrational excitations in the molec- 
ular region. Also there were some efforts to include dynamical correlations 
into DFT-based approach. Much harder and till the present not achieved at all 
is the inclusion of electronic correlation effects, responsible for many-particle 
effects like Coulomb blockade or the Kondo effect, in a non-equilibrium trans- 
port situation. DFT-based techniques, being essentially mean-field theories, 
cannot deal in a straightforward way with such problems and have to be 
improved, e.g. within the LDA+U approaches. For the case of equilibrium 
transport, a generalization of the Landauer formula including correlations 
has been recently formulated as well as first attempts to go beyond the lin- 
ear response regime; for strong out-of-equilibrium situations this will be, in 
our view, one of the most demanding issues that the theoretical "transport" 
community will be facing in the coming years. 
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